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TECHNICAL  REPORT  SUMMARY 


This  report  covers  contract  work  for  the  time  period  01  Oct  1978  to  30  Apr 

ductlve  that  I  have  spent  on 
he  report  Is  divided  Into 
two  parts.  The  first  part  describes  contract  work  that  Is  now  being  written  up 
In  publishable  manuscripts.  Topics  covered  Include:  Ill  analysis  of  wideband 
digital  data  for  elucidation  of  the  seismic  source  (two  manuscripts  In  prepara¬ 
tion),  121  analysis  of  eastern  Kazakh  explosions  for  the  effects  of  anomalous 
tectonic  excitation  (1  In  preparation),  and  13)  an  evaluation  of  the  seismic  mo¬ 
ment  tensor  method  for  analysis  of  close-in  accelerometer  data  of  the  JORUM, 
PIPKIN,,  and  HANDLEY  tests  (1  in  preparatl on).  The  second  part  of  the  report 
consists  of  4  manuscripts  that  have  been  completed  for  submission  to  the  jour¬ 
nals.  Topics  covered  Include:  14)  amplitude-yield  scaling  of  chemical  explo¬ 
sions;  15)  new  evidence  on  the  Pn-Raylelgh  and  another  neat — regional  seismic 
discriminant,  16)  the  computation  of  near-field  canonical  sources  by  exact  theo¬ 
ry  in  a  multilayered  medium,  and  17)  seismic  discrimination  at  near-regional 
distances  using  only  Pn. 


198X* ■  The  time  has  been  the  most  successful  and 
DARPA/AFOSR  objectives  in  seismic  discr Imlnatlon 
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(1)  ANALYSIS  OF  WIDEBAND  DIGITAL  DATA 


Ever  since  my  participation  In  the  Near-Field  Project  sponsored  by  AFOSR,  I 
have  felt  that  a  primary,  ongoing  need  for  the  advancement  of  earthquake  source 
theory  is  the  acquisition  of  better  data  than  had  hitherto  been  available; 
source  models  are  numerous,  but  seismologists  still  disagree  on  some  of  the  most 
fundamental  points.  Therefore,  I  have  felt  a  responsibility  to  obtain  such  data 
as  time  and  resources  have  permitted.  With  the  advent  of  digital 
event-recording  seismographs,  ve  began  a  program  of  digital  data  acquisition  on 
eastern  Sierra  earthquake  sequences.  The  rate  of  data  return  has  been  very  high 
compared  with  the  modest  costs  of  acquiring  it. 


We  were  fortunate  in  that  eastern  Sierra  seismicity  began  to  Increase,  af- 
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fording  the  opportunity  to  record  a  number  of  distinct  aftershock  sequences.  We 
were  extremely  fortunate  to  have  placed  event  recorders  In  the  field  6  months 
prior  to  the  onset  of  the  major  sequence  at  Mammoth  Lakes,  culminating  In  ten  ML 
5+  events  in  the  last  week  of  May,  1980.:  The  area  Is  of  considerable  tectonic 
interest,  and  these  sequences  led  us  to  a  number  of  research  projects  (Somei — 
ville  and  Peppin,  1980a, b;  Peppin  and  Somerville,  1980;  Peppin  and  Somerville, 
1981a).  The  work  relevant  to  DARPA  was  done  on  an  extension  of  the  original 
contract,  and  Involved  processing  the  digital  data  acquired.  This  work  will  be 
described  in  two  manuscripts  In  preparation  (Peppin  and  Somerville,  1981b; 
Somerville  and  Peppin,  1981). 

A.  The  digital  data  set 

In  Figure  is  presented  an  image  showing  the  locations  and  number  of 
three-component  digital  seismograms  recorded.  The  data  is  recorded  wideband 
displacement  flat  (0.1  to  50  Hz)  using  seismograph  systems  as  described  by  Pep¬ 
pin  and  Bufe  (1980).  Altogether  some  1,500  3-component  seismograms  were  ana¬ 
lyzed.  For  almost  all  events  analysis  consisted  of:  (1)  dumping  the  digital 
data  to  named  mass-storage  files;  (2)  computing  amplitude  spectra  of  the  verti¬ 
cal-component  P  and  horizontal-component  S-waves;  (3)  constructing  plots  of  all 
data  traces  and  spectra;  (4)  measuring  spectral  corner  frequencies;  (5)  com¬ 
puting  seismic  moment  from  the  long-period  spectral  asymptote;  and  (6)  comput¬ 
ing  synthetic  Wood-Anderson  magnitudes.  In  addition  we  have  computed  adjusted 
Berkeley  magnitudes  where  possible  and  seismic  stress  drop  using  Brune's  (1970, 
1971)  formulas.  More  detailed  analysis  awaits  obtaining  accurate  locations  for 
study  events  of  interest,  a  problem  rendered  nontrivial  for  several  reasons, 
mostly  because  of  complex  geology,  poor  station  distribution,  and  record  overlap 
in  foreshocks.  Data  quality  is  superior,  with  large  s igna 1-to-no 1 se  ratio  ex¬ 
tending  over  most  or  all  of  the  recording  bandwidth. 

B.  Pictorial  description  of  the  data  set 
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For  so  large  a  data  set.  It  Is  of  Interest  to  provide  some  sort  of  visual 
characterization.  Presented  here  are  plots  of  log  seismic  moment  versus  ML 
(F Iqures  1.2  and  1.3),  log  of  seismic  moment  versus  P-  and  S-wave  spectral 
corner  frequency  (F loures  1 . 4  and  1.6).  log  of  seismic  moment  versus  their  ratio 
(F i gure  1.6).  and  their  ratio  versus  P-wave  spectral  corner  frequency  (Figure 
1.7). 


C.  Discussion 

The  original  purpose  of  this  experiment  was  to  acquire  a  set  of  excellent 
data  close  to  the  source  of  earthquakes.  Then  we  could  experiment  with  waveform 
Inversion  methods  for  determining  the  seismic  source  function.  One  problem  In 
particular  remains  in  question  today,  namely!  what  source  parameter  relates  to 
the  time  duration  of  seismic  waves  eminating  from  the  source  (or  similarly,  to 
the  spectral  corner  frequencies  of  P  and  S  waves)?  The  prevailing  view  has  been 
that  the  spatial  source  dimension  (i.e.,  length  of  faulting)  controls  the  ob¬ 
served  frequencies.  However,  Helmberger  and  colleagues  have  had  good  success 
modelling  moderate  earthquakes  essentially  as  point  sources  in  their  studies  of 
teleselsmic  body  waves.  Indeed,  for  our  data  set  it  is  remarkable  how  far  the 
point-source  approximation  can  be  taken.  In  F Iqure  1 .8  we  show  data  recorded 
close  (1  second  6  -  P  time)  to  two  sizeable  Mammoth  aftershocks.  Synthetic  se¬ 
ismograms  for  a  point  source  in  a  halfspace  have  also  been  prepared,  and  the 
fits  are  remarkable  under  the  circumstances.  The  larger  earthquake  has  a  seism-- 
1c  moment  of  1.7E23  dyne-cm,  one  of  the  bigger  events  In  the  sequence.  And  yet, 
as  can  be  seen,  the  effect  of  source  finiteness  Is  difficult  to  identify. 

Hanks  (1981)  explains  how  the  point-source  assumption,  If  Incorrect,  could 
lead  to  distinctly  wrong  values  of  mantle  Q  for  P  (versus  S)  waves.  The  reason 
for  the  problem  is  this.  What  if  In  fact  it  is  a  general  result  that  the 
P-waves  leaving  the  source  are  richer  In  high  frequencies  than  the  S-waves. 
Then  we  would  expect  to  see  this  effect  in  teleseismlc  body  waves  (which  we  do). 
But  if  the  difference  is  then  attributed  to  whole-path  Q,  we  will  erroneously 
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Figure  1.4.  Log  seismic  moment  versus  log  of  the  vertical -component 
P-wave  spectral  corner  frequency.  These  were  inferred  by  overlaying 
templates  on  the  spectra  and  picking  the  corners  by  eye,  most  to  with¬ 
in  a  confidence  range  of  about  20%. 
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Figure  1 .5.  Log  of  seismic  moment  versus  log  of  horizontal -component 
S-wave  spectral  corner  frequency.  Inferred  by  overlaying  a  template 
and  picking  by  eye.  Most  are  precise  to  about  20%. 
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Figure  1.6.  Log  of  seismic  moment  versus  the  ratio  "R"  of  vertical- 
component  P  versus  horizontal -component  S  spectral  corner  frequency 
Because  of  overlapping  points,  it  is  not  clear  in  the  figure  that  the 
average  of  all  values  is  about  1.3. 
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Fi gure  1.7.  Spectral  corner-frequency  ratio  "R"  versus  vertical -component 
P-wave  spectral  corner  frequency  at  all  stations  except  MCE  and  MRS  (top), 
MRS  (lower  left),  and  MCE  (lower  right).  Dashed  lines  are  trend  lines  pre¬ 
dicted  as  described  in  text.  Numerals  beside  each  line  indicate  cutoff 
frequency  for  S  selected.  At  all  stations  except  MRS  there  seems  to  be  a 
noticable  correlation  betv/een  R  and  the  P-wave  spectral  corner  frequency. 
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Figure  1 .8.  Observed  and  theoretical  seismograms  for  two  large  events,  26  May  1925, 

ML  4.6  (left,  transverse)  and  26  May  1980  1857,  ML  5.5  (vertical  and  transverse, 
right).  Ramps  between  P  and  S  phase  are  the  near-field  terms,  visible  on  vertical 
as  well  as  transverse  components.  For  1857  the  preamplifiers  are  clipped  after  onset 
of  the  S-phase  (last,  large  arrival  on  all  traces). 
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Infer  values  of  Q  for  P  too  high  relative  to  Q  for  S. 

Now,  says  Hanks,  we  do  see  this  frequency  shift  In  the  majority  of  data 
studies  published,  regardless  of  recording  distance.  Therefore,  he  asserts  that 
the  shift  is  a  source  effect,  which  implies  that  the  length  of  the  rupture  zone 
controls  the  frequencies  observed. 

With  the  sizeable  data  set  at  hand,  and  given  its  superb  quality,  Important 
new  evidence  can  be  brought  to  bear  on  this  question.  In  Figure  1.6  we  have 
plotted  the  ratio  "R"  of  P  to  S  wave  spectral  corner  frequency  for  the  events 
possessing  good-quality  corner  frequencies  on  each  phase.  We  see  that  the  data 
support  Hanks  in  that  the  P-waves  appear  to  be  enriched  in  h i gh-frequency  energy 
relative  to  S.  However,  we  can  give  a  different  1 nterpretat  . on .  At  each  re¬ 
cording  station  a  "cutoff  frequency"  can  almost  always  be  identified  for  the 
S-waves  spectral  corners,  but  not  the  P.  This  cutoff  is  frequently  quite  clear¬ 
ly  identifiable.  Consider,  for  example,  the  data  at  the  Whitmore  station  shown 
in  F i qure  1.9.  The  S-wave  data  seem  cut  off  at  all  levels  of  seismic  moment  at 
about  7  Hz,  while  P-wave  corners  go  up  to  15  Hz  or  more  for  fairly  large  events. 
One  obvious  explanation  for  these  facts  is  the  action  of  anelastic  attenuation. 
Suppose  that  this  effect  works  in  such  a  way  as  to  band-limit  S-waves  while 
leaving  P  unaffected.  Then  we  would  expect  R  to  correlate  with  P-wave  spectral 
corner  frequency. 

In  Figure  1.7  we  suppose  that  R  is  unity  at  the  source  and  that  a  separate 
cutoff  frequency  exists  at  each  station  for  S.  Above  the  cutoff  the  observa¬ 
tions  of  R  should  begin  to  increase.  The  lines  in  Figure  1.7  are  predicted 
trend  lines  on  the  assumption  of  varying  cutoff  frequencies  for  S;  the  fit  to 
the  observations  is  reasonable  for  cutoff  values  close  to  those  as  observed  in 
Figure  1.9. 


From  this  analysis  we  conclude  that  no  statement  can  be  made  about  condi¬ 
tions  at  the  source  based  on  the  observations  (Figures  1.4  and  1.7)  until  the 
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Figure  1.9.  S-wave  spectral  corner  frequency  (abscissa)  versus 
log  of  seismic  moment  for  events  recorded  at  the  Whitmore  station. 
The  corner  frequency  data  are  almost  cutoff  completely  at  7  Hz. 
Such  cutoffs  can  be  identified  at  all  of  the  stations  except  DIA. 
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effects  of  anelastic  attenuation  are  better  known,  so  that  Hanks’  (1981)  premise 
Is  not  justified  for  this  data  set.  In  fact,  let  cutoff  frequencies  be  selected 
from  Figure  1.7.  Then  let  R  be  computed  only  when  the  P-wave  corner  Is  strictly 
not  greater  than  this  cutoff.  At  8  of  9  stations,  the  average  value  of  R  is  in¬ 
distinguishable  from  unity  (Table  1).  In  support  of  this  general  idea,  note 
that  the  rate  of  falloff  to  high  frequency  is  observed  to  be  significantly  gre¬ 
ater  in  S-waves  than  in  P-waves,  consistent  with  the  proposition  that  the 

S-waves  are  being  more  strongly  affected  by  anelastic  loss.  In  contrast,  I  know 
of  no  source  theory  that  can  predict  such  a  difference  in  high-frequency  falloff 
on  elastic  rheology. 

To  me  all  this  means  that  R  is  probably  unity  at  the  source  for  these 

events.  Then  it  must  be  the  source  rise  time  that  dictates  the  spectral  corner 

frequencies,  not  the  spatial  source  dimension.  If  true,  this  result  can  also 

explain  the  observation  of  rather  weak  scaling  of  the  P-wave  spectral  corner 
frequencies  (factor  of  near  2  over  4  orders  of  magnitude  in  seismic  moment). 
Source  dimensions  should  scale  more  rapidly  with  moment  than  this  on  any  reason¬ 
able  hypothesis. 

I  believe  these  results  are  quite  relevant  to  the  DARPA  discrimination  pro¬ 
gram,  because  they  relate  in  a  fundamental  way  to  our  treatment  of  earthquake 
sources.  I  can  see  no  way  to  make  sense  of  the  observations  unless  these  earth¬ 
quakes  are  strongly  controlled  at  all  magnitude  levels  by  asperities  of 
more-or-less  limited  size  range  (say  25  ±  10  m>.  In  this  case,  high-frequency 
body  waves  are  seeing  the  failure  of  rather  small  asperities,  while  surface 
waves  are  seeing  the  effects  of  a  rupture  that  radiates  out  over  a  larger,  but 
less  energy-rich  rupture  surface.  "Stress  drop"  for  such  a  failure  has  no  use¬ 
ful  meaning  when  taken  over  the  whole  fault,  but  might  have  meaning  correlable 
with  the  strength  and/or  texture  of  rocks  in  the  source  region.  If  this  Is  a 
general  result,  it  means  that  theoretical  treatment  of  earthquake  sources  for 
the  discrimination  problem  should  be  modified  in  fundamental  ways;  for  earth¬ 
quakes  are  then  highly  "dissimilar”  in  the  sense  of  Aki  (1967).  If  so,  our  dis- 
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Table  1.  Spectral  Corner  Frequency  Ratios 


Station 


MCE 

MRS 

ROC 

WIT 

BIR,BIW 

BRB.BRS 

DIA 


Sequence  R(all  data)  N  cutoff  freq.  R  (P  cutoff)  N 


Mammoth  1.17  t 
Mammoth  1.18  - 
Mammoth  1.33  ± 
Mammoth  1.61  * 
Bishop  1.81  t 
Bridgep.  1 .62  ± 
Dia.  Val .  2.30  ± 


.29  91  13 

.28  152  8 

.31  101  8 

.53  160  7 

.66  37  8 

.58  52  8 

.40  18 


0.98  +  .26  30 

1.15  +  .33  81 

1.16  ±  .26  44 

1.26  ±  .27  62 

1.12  t  .18  13 

0.96  ±  .09  11 
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crimination  efforts  should  be  directed  more  toward  statistical  measures  of  this 
dissimilarity  in  source  regions  of  Interest  rather  than  the  present  heavy  em¬ 
phasis  on  waveform  fitting  using  relatively  simple  sources.  Where  dlsslml lar Ity 
is  great,  scatter  (hence  overlap)  of,  say,  Ms:mb  populations  of  explosions  and 
earthquakes  can  be  expected. 

C2)  SEISMIC  CHARACTERIZATION  OF  RUSSIAN  EVENTS 

Of  recent  Interest  to  the  DARPA  program  was  the  occurrence  of  the  event  of 
07  July  1979,  an  eastern  Kazakh  presumed  explosion  showing  a  tremendous  SH  pulse 
at  the  SRO  station  ANTO  (At  a  recent  DARPA  contractors'  meeting,  T.  Fitch  ex¬ 
pressed  the  opinion  that  this  pulse  was  not  real;  If  so,  then  It  Is  notorious 
that  the  motion  is  almost  pure  transverse  on  the  two  horizontal  components  rela¬ 
tive  to  the  great-circle  travel  path).  Brian  Stump  and  I  were  asked  to  analyze 
this  event  with  a  view  toward  answering  the  following  3  questions:  (a)  can  the 
event  he  identified  as  an  explosion?  (b)  can  the  tectonic  orientation  of  the 
stresses  involved  in  production  of  the  SH  phase  be  identified?  (c)  can  we  esti¬ 
mate  the  effects  of  the  added  earthquake  component  on  the  amplitude  of 
short-period,  vertical  P  (i.e.,  the  effect  on  the  body  wave  magnitude)? 

Invest igat i ons  completed  so  far  can  address  (a)  and  (c)  and  Indicate  ways 
to  proceed  toward  (b).  We  began  by  establishing  a  control  event,  that  o^  08  Au¬ 
gust  1979,  for  which  the  best  short-period  data  was  available,  and  for  which  no 
large  SH  pulse  was  evident  on  the  SRO  stations.  Using  eight  SRO  and  ASRO  sta¬ 
tions,  we  inverted  for  the  seismic  moment  tensor.  For  a  Eurasian  mantle  struc¬ 
ture  provided  by  Helmberger  and  WKBJ  analysis  (Chapman,  1373;  Dey~Sarker  and 
Chapman,  1978)  for  the  Green's  functions,  the  following  extreme  (unsigned)  va¬ 
lues  of  the  moment  tensor  Mij(t)  were  obtained: 
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A  pure  explosion  would  have  equal  values  on  the  diagonal  and  zeroes  off  the 
diagonal.  Thus,  the  event  is  inferred  to  be  part  explosion  and  part  earthquake. 
We  consider  the  result  a  success  because:  (1>  only  vertical  component  P  data 
was  usedj  (2)  good  fits  of  the  data  at  all  stations  and  somewhat  reliable  time 
functions  were  obtained;  (3)  the  phase  sP  was  not  included  in  the  Green’s  func¬ 
tion  due  to  a  bug  in  the  WKBJ  code  (inclusion  of  this  phase  would  no  doubt  drive 
down  the  off-diagonal  elements  that  are  very  sensitive  to  S  leaving  the  source). 
For  samples  of  the  fits  of  the  data,  see  Figure  2.1,  and  for  samples  of  the  in¬ 
ferred  motions  at  the  source,  see  F i gure  2.2. 

Regarding  the  seismic  discrimination  question  (a):  it  appears  likely  that 
the  gross  nature  of  the  source — explosion  or  earthquake — can  be  resolved  by  this 
procedure  using  teleseismic  body  waves,  although  investigations  of  other  explo¬ 
sions  and  earthquakes  will  be  required  to  verify  this.  Our  conclusion  regarding 
(c)  is  that  mb,  hence  inferred  yield,  is  unlikely  to  be  affected  by  "anomalous’ 
tectonic  excitation  at  the  level  of  the  07  July  1979  event.  Question  (b)  rema¬ 
ins.  Recent  communications  with  C.  H.  Chapman  have  verified  some  mistakes  In 
the  WKBJ  source  terms  employed  by  Stump  and  me  for  eastern  Kazakh  events.  Thus, 
we  need  to  recompute  the  necessary  Green's  functions  incorporat ing  the  phase  sP. 
For  the  08  August  1979  event  we  were  able  to  infer  principal  directions  of 
stress  for  the  deviatoric  part  of  the  solution.  The  results  was  a  thrust  me¬ 
chanism  with  compression  axis  at  N  15  deg  E,  plunge  10  degrees,  and  with  tension 
axis  at  about  S  20  deg  E,  plunge  50  degrees.  These  Inferred  directions  are 
likely  to  be  altered  signif Icantly  by  inclusion  of  sP  in  the  Green's  functions, 
but  our  other  conclusions  about  these  sources  will  probably  not  be  altered  sub- 
stant i ve ly . 

131  CLOSE-IN  ACCELEROMETER  DATA 

The  University  of  California  at  Berkeley  conducted  experiments  to  record 
ground  motions  near  several  large  underground  explosions  JORUM  (16  Sept  1969), 
PIPKIN  (08  Oct  1959),  and  HANDLEY  (26  Mar  1970).  Results  of  the  analyses  of 


Figure  2.1.  Data  (top  trace)  and  fit  (bottom  trace  of  each  pair)  that 
results  from  fitting  8  vertical -component  short-period  SRO  and  ASRO 
records  of  the  eastern  Kazakh  explosion  of  08  Aug  1979.  Fits  for  stations 
not  shown  are  as  good,  except  the  CHTO  record.  Distances  in  degrees  ("d") 
and  station-to-epicenter  azimuth  ("az”)  are  given.  Traces  are  scaled  for 
relative  amplitude  between  data  and  fits. 
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Figure  2.2.  Inferred  motions  at  the  source  for  the  eastern  Kazakh  explosion 
shown  in  Figure  2.1.  Top:  trace  of  the  moment  tensor  (explosive  part); 
Middle:  M-j^t);  Bottom:  MopCt).  Note  the  prominent  arrival  after  the 
onset  suggestive  of  spall  closure.  Note  the  roughly  steplike  appearance 
of  the  deviatoric  component  (center  trace).  Finally,  note  that  these  records 
(top  trace)  seem  to  indicate  that  the  explosive  part  of  the  source  has  over¬ 
shoot. 
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these  records  are  a  major  part  of  two  theses  (Peppln,  1974;  Stump,  1979),  as 
well  as  recent  work  by  Helmberger  and  Hadley  (1981)  and  Stump  and  Johnson 
(1981).  The  aim  here  has  been:  determine  the  source  function  of  nuclear  explo¬ 
sions  using  observations  taken  close  to  the  source;  then  questions  pertinent  to 
seismic  discrimination,  e.g.,  the  shape  of  the  time  function  and  amplitude-yield 
scaling,  will  improve  our  understanding  of  why  the  various  discriminants  are  ef¬ 
fective  . 

Recently  Brian  Stump  and  I  put  our  two  codes  together.  My  code  produces 
Green's  functions  for  multi  layered  media  close  to  the  source  with  no  approxima¬ 
tion,  having  been  designed  specifically  for  this  problem  (Peppin,  1981),  this 
report),  and  his  does  the  linear  inversion  for  the  source.  I  have  made  four 
runs  on  the  HANDLEY  data  and  one  with  the  PIPKIN  data  with  the  aim  of  addressing 
these  questions:  how  stable  are  the  results  to  large-scale  perturbations  in 
structure  and  what  can  we  infer  as  true  source  motion? 

The  answer  is  mixed.  Stump  and  Johnson  (1981)  have  considered  the  effect 
of  computing  Green's  functions  for  one  structure  with  different  sets  of  general¬ 
ized  rays,  while  here  we  consider  the  effects  of  distinctly  different  layer  mo¬ 
dels.  If  we  make  a  comparison  based  on  the  maximum  values  of  the  moment  tensor, 
as  presented  in  the  previous  section,  then  all  structures  shown  In  F 1 qure  3. 1 
lead  to  the  unequivocal  character i zat 1  on  of  the  event  as  dominantly  an  explosion 
as  found  by  Stump  and  Johnson.  The  models  for  the  halfspace  and  for  the  source 
in  the  top  layer  (Figure  3.1)  show  a  source  which  exhibits  departure  from  sphei — 
ical  symmetry  in  the  vertical  direction;  the  Helmbergei — Hadley  and 
source-below-tne-layer  models  in  Figure  3.1  give  values  of  M(1,J)  with  the 
off-diagonal  terms  a  twentieth  or  less  of  the  on-diagonal  terms  and  these  are 
within  a  few  percent  of  being  identical  in  magnitude  and  shape.  All  of  the  mo¬ 
dels  in  Figure  3.1  lead  to  excellent  fits  of  all  12  components  of  recorded  ac¬ 
celeration  on  vertical,  radial,  and  transverse  components  essentially  throughout 
the  whole  record.  As  a  result  of  these  efforts  I  am  convinced  that  one  signifi¬ 
cant  problem  has  been  identified  and  solved.  For  the  halfspace  Green's  function 
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Figure  3.1 .  Layer  models  used  in  source  inversions  of  close-in  accelerometer  data. 
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or  for  the  Green's  function  corresponding  to  a  source  In  the  top  layer,  the  the¬ 
oretical  solutions  predict  significantly  larger  radial  than  vertical  motion 
while  the  reverse  is  observed  in  the  main  P-wave  of  the  explosion.  This  leads 
directly  to  the  inference  of  an  axial  component  in  the  source  (the  inversion  has 
no  other  way  to  make  the  vertical  component  larger  than  the  radial).  Thus,  this 
axial  component  is  definitely  an  artifact,  not  part  of  the  source.  When  any 
Green's  function  is  used  which  gives  the  correct  vert i ca 1 /rad  1  a  1  amplitude 
ratio,  the  axial  component  disappears,  leaving  a  spher i ca  1  ly-symmetr i c  source. 
Consequently,  I  am  now  inclined  to  reject  the  non-explosive  source  component  In¬ 
ferred  by  Peppin  (1977),  at  least  insofar  as  the  generation  of  close-in  body 
waves  is  concerned. 

The  problem  is  that  the  inferred  time  histories  of  the  source  are  quite 
variable  from  model  to  model;  compare  the  time  histories  inferred  for  the  ex¬ 
plosion  part  of  the  source  presented  in  F iqure  3.2.  It  appears  that  the  method 
is  such  that  the  relative  strengths  of  the  M<i,j>  are  insensitive  to  changes  in 
the  model,  but  the  phase  information,  which  controls  the  shape  of  the  time  his¬ 
tory,  is  sensitive.  That  is  to  say:  we  can  get  some  useful  information  with 
structure  known  only  vaguely,  but  we  can  get  much  more  if  the  structure  is  known 
to  better  precision.  Stump  and  Johnson  (1981)  find  elements  in  common  with  all 
of  their  inferred  sources;  however,  their  solutions  appear  unsatisfactory  to  me 
in  that  the  most  energetic  part  of  the  source  lags  1.3  seconds  behind  the  first 
impulse.  More  work  is  needed  on  this  problem. 

Cnee  satisfactory  inversions  can  bo  perforn?d,  our  aim  is  to  try  and  find 

elements  of  similarity  in  the  source  functions  inferred  under  different  circum¬ 
stances  with  different  explosions  (do  all  need  overshoot?  do  all  need  no  ovei — 
shoot?).  Vie  want  to  know  whether  a  simple  source  is  sufficient,  or  if  spall 
phenomena  contribute  (here  note  the  second  arrival  In  the  source  function  about 
2  seconds  after  the  event  onset  in  Figure  3.2,  correspond i ng  to  spall  closure 
time  for  an  event  of  this  magnitude:  Springer,  1974).  The  overshoot  question 
has  been  studied  using  the  same  data  set  by  Peppin  (1977),  Stump  (1979),  Helm- 
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Figure  3.2.  Inferred  source  motions  for  the  explosive  part  of  the  source 
of  the  HANDLEY  explosion.  Top  trace:  for  the  model  left-center  in  Figure 
3.1;  bottom  trace:  for  the  model  right-center  in  Figure  3.1.  Prominent 
secondary  arrival  seen  on  second  trace  as  found  by  Stump  (1979)  for  a  half¬ 
space  and  by  Stump  and  Johnson  (1981)  for  the  bottom  layer  model  of  Figure 
3.1.  Again  the  source  appears  to  have  overshoot;  compare  with  Figure  2.2. 
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berger  and  Hadley  (1981)>  and  Stump  and  Johnson  (1981)  with  diverging  opinions. 
We  must  be  prepared  to  accept  that  different  source  functions  may  be  required  to 
explain  data  at  close-in  versus  teleseismic  distances;  and  we  must  then  assess 
the  implications  for  seismic  discrimination  insofar  as  transferring  Inferences 
based  on  the  one  data  set  to  the  other. 
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[41  AMPLITUDE-YIELD  SCALING  RELATIONS  OF 
CHEMICAL  EXPLOSIONS 


William  A.  Peppin 


Vaughn  E.  Culler 


SUMMARY 

Forty-six  chemical  explosions  (24  to  515  pounds  of  C4  equivalent)  were  de¬ 
tonated  above  ground  in  a  caref u 1 ly-contr o 1  led  experiment  for  seismic 
anpl i  iude-y is  Id  scaling  relations.  Wideband  (0.1  tc  50  Hz)  digital  seismic  data 
was  recorded  at  four  sites  1-4  km  from  various  shotpolnts.  The  seismic  phases 
analyzed  were  a  P-wave  refracting  along  a  shallow  interface  and  the  ground  res¬ 
ponse  to  the  air  blast.  Twenty-seven  amplitude-yield  scaling  parameters  were 
developed.  Analysis  of  analog  and  spectral  parameters  shows  the  superiority  of 
the  latter  in  several  respects:  reduced  dependence  on  source  coupling  and  re¬ 
ceiver  conditions;  greater  case  of  measurement  standardization;  improved  ob¬ 
jectivity  in  making  a  measurement.  These  results  encourage  the  investigation  of 
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spectral  amplitude-yield  scaling  relations  for  the  phase  Pn  at  regional  distance 
ranges . 


INTRODUCTION 

A  sizeable,  literature  exists  on  the  subject  of  seismic  amplitude-yield 
scaling  (Carder  and  Cloud,  1 959 ;  Romney,  1959;  Adams  et  al.,  1961;  Carder  and 
Mickey,  1962;  Kovach  et  al.,  1963;  Murphy  and  Lahoud,  1969;  Springer  and  Han¬ 
non,  1973).  However,  summarized  results  (Murphy  and  Lahoud,  1969,  Tables  2  and 
3;  Springer  and  Hannon,  1973)  differ  markedly  with  charge  size,  recording  In- 
strumentat  i  on  ,  distance  range,  and  source/receiver  geology.  At  Lawrence  Llvei — 
more  Laboratory  (LLL),  it  was  felt  that  insufficient  data  was  available  In  the 
relevant  ranges  of  charge  size  and  recording  distance.  Accordingly,  we  designed 
a  seismic  experiment  which  was  conducted  at  the  Site  300  testing  facility  near 
LLL  for  the  purpose  of  developing  a  set  of  ampl Itude-yield  scaling  relations  for 
the  seismic  amplitudes  produced  by  above-ground  chemical  explosions. 

Two  facts  establish  the  distinctive  nature  of  our  experiment  as  compared 
with  previously-conducted  ones.  First,  we  collected  some  of  the  best  seismic 
data  ever  taken  for  this  purpose,  and  second,  the  design  of  the  experiment  was 
left  to  a  seismologist  within  fairly  broad  constraints.  Thus,  we  were  able  to 
obtain  exactly  the  data  we  thought  relevant,  and  because  this  data  was  wideband, 
three-component,  velocity-flat  (0.1  to  50  Hz;  12-bit  resolution;  X16 
gain-range  step),  complete  waveforms  were  available  for  all  charges  and  record¬ 
ing  distance  ranges.  A  detailed  description  of  the  seismograph  systems  employed 
(Dr.Rs)  is  found  in  Puppin  and  Bufe  (1980). 

EXPERIMENTAL  PROCEDURE  AND  RECORDS  TAKEN 

Three  DERs  were  deployed  1  km  SSE,  3  km  SE ,  and  2.3  km  SSW  of  the  standard 
802  shotpoint.  These  were  given  respective  station  names  845,  Llnac,  and  858 


(Figure  1).  Threc-c omponent  data  was  recorded  at  each  site  along  with  WWVB,  an 
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absolute  time  standard.  Horizontal  components  were  rotated  to  radial  (positive 
away)  and  transverse  (positive  to  the  right  facing  the  source)  components. 
Tests  were  made  which  showed  Independence  of  amplitude  scaling  with  charge  he¬ 
ight  above  the  ground  (0  to  2  m)  and  with  respect  to  asymmetric  source  modifica¬ 
tion  (200  kg  steel  plates  placed  above,  beside,  and  below  100-pound  charges). 
For  the  amplitude-  yield  scaling  series  “KSG",  each  charge  was  fired  30  cm  above 
the  ground,  the  crater  was  refilled,  and  the  next  shot  fired  In  the  same  spot  to 
a  precision  of  4  m.  Yields  ranged  from  24  to  515  pounds  of  C4  equivalent  high 
explosive.  In  the  source  coupling  series  "KSF11  Identical  50-pound  charges  were 
fired  In  an  X  pattern  over  the  802  shotpoint  with  KSF4  at  that  site  for  refei — 
ence  (Figure  1). 

Figure  2  shows  a  typical  record  as  seen  at  the  Linac  site.  On  these  and 
all  other  traces  two  phases  are  seen.  The  first  motion  is  a  seismic  P-wave 
(longitudinal  from  the  source,  hereinafter  called  the  "seismic  wave")  and  the 
second  is  the  ground  response  to  the  air  blast  (hereinafter  the  "acoustic 
wave"),  which  was  fairly  loud  and  sensible  at  this  site.  Figure  3  is  a 
travel-time  plot  for  the  seismic  wave.  The  velocity  Is  3.06  km/s  and  the  inter¬ 
cept  is  consistent  with  a  shallow  refractor  (typical  at  Site  300  is  60  m  of 
loose  material  over  Cenozoic  sedimentary  rocks).  Thus,  the  arrival  undoubtedly 
is  a  refraction  along  the  top  of  the  Cenozoic  rocks,  analogous  to  Pn  observed 
regionally  in  the  western  U.S.  In  Figure  2  note  the  long-period  undulations 
seen  on  each  component  preceding  the  event  onset.  These  are  the  six-second  mi¬ 
croseisms,  which  demonstrates  the  system  recording  bandwidth.  We  identified  no 
other  seismic  phase  (such  as  S,  Love,  or  Rayleigh  waves).  Peppin  (1978)  givos  a 
complete  presentation  of  all  data  recorded  and  a  fuller  discussion  of  these  re¬ 
cords  . 


DATA  ANALYSIS 


Analog  measurements  were  made  from  the  400  traces  available  from  each 

component  of  seismic  motion  recorded.  Amplitude  spectra  of  the  seismic  and 
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acoustic  phase  (800  In  all)  were  computed.  The  windows  were  1.28  and  2.56  sec¬ 
onds,  respectively,  for  these  phases.  A  least-squares  parabolic  trend  was  re¬ 
moved  and  a  10%  cosine  taper  applied  at  each  end.  Spectra  were  computed  from 
the  resulting  traces.  These  were  smoothed  (3-polnt  triangle  window)  and  all 
were  plotted  after  converting  to  ground  motion  by  removing  the  system  response. 
Identical  processing  was  applied  to  the  quiet  segment  preceding  event  onset  so 
that  an  estimate  of  slgnal-to-noise  ratio  could  be  obtained  as  a  guide  for  se¬ 
lecting  appropriate  passbands  for  amplitude-yield  scaling.  A  sample  spectrum 
appears  in  Figure  4;  the  others  are  found  in  Peppin  (1978).  All  of  the  scaling 
measurements  were  made  either  from  the  preprocessed  time  trace  (lower  part  of 
Figure  4)  or  from  a  spectral  average  of  selected  passbands  In  the  amplitude 
spectrum  of  this  trace  (upper  part  of  Figure  4).  The  analog  measurements  taken 
were  the  first  quarter-cyc le  amplitude  ("a")  and  first  half-cycle  amplitude 
("b" )  measured  directly  off  the  time  trace  in  Figure  4  (all  are  plotted  the  same 
height  and  the  maximum  amplitude  in  digital  counts  is  given  to  the  right  of  the 
time  trace).  Spectral  scaling  averages  were  determined  by  trial  and  error  (at 
Linac  the  band  6-9  Hz  was  selected).  "Better"  passbands  were  those  giving  more 
uniform  amplitude-yield  scaling  of  the  KSG  series.  We  claim  no  completeness  In 
our  search  for  possible  scaling  parameters. 

RESULTS 


A.  Source/Recei ver  Effects 

The  shots  of  the  KSG  series,  all  on  802,  provided  seismic  records  which 
overlaid  site  by  site  and  component  by  component  apart  from  a  scaling  of  ampli¬ 
tudes  with  charge  size.  Thus,  amplitude-yield  scaling  was  precise  to  5%  and  ac¬ 
curate  to  10%  for  charges  fired  at  that  site  and  recorded  at  any  of  the  three 
seismographs.  The  KSF  series  was  designed  to  investigate  variable  source  cou¬ 
pling  (see  Figure  1).  In  Table  1  we  consider  the  ratio  of  a  given  measurement 
resulting  from  a  shot  off  802  to  the  same  measurement  resulting  from  the  KSF4 
shot  on  802.  Since  all  charges  of  this  series  had  the  same  weight  (50  pounds), 
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Table  1.  Variability  in  source  coupling  by  measurement  type 


Quantity  Measured 


Fraction  of  KSX4 


linac  seismic,  spectral 
845  seismic,  spectral 
858  seismic,  spectral 


.753 

.944 

.609 


Linac  acoustic,  spectral 
845  acoustic,  spectral 
858  acoustic,  spectral 


Linac  seismic,  analog 
845  seismic,  analog 
858  seismic,  analog 


.652 

.635 

.762 
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a  unit  ratio  Implies  Identical  source  coupling.  From  Table  1  note  the  signifi¬ 
cant  variation  In  source  coupling  (shots  spaced  only  tens  of  m  apart)  and  that 
spectral  measurements  are  less  source-dependent  than  analog  ones.  Note  also 
that  the  acoustic  wave  shows  no  variation  In  source  coupling,  as  Is  expected. 

Table  2  provides  a  more  detailed  breakdown  for  each  shot  of  the  KSF  series. 
In  this  table  we  present  the  ratios  for  the  seismic  analog  and  spectral  measure¬ 
ments  by  shot.  In  column  4  we  give  the  ratio  we  would  expect  and  In  column  5  we 
state  why;  for  example,  as  KSF2  was  placed  on  a  dry  sandstone  outcrop,  we  would 
expect  reduced  motions  near  the  source  as  compared  with  KSF4  fired  In 
loosely-consolidated  mud  and  gravel.  This  table  shows  that  in  each  case  except 
KSF8  the  spectral  measurement  gave  the  more  reasonable  ratio  relative  to  KSF4. 

To  pursue  these  effects  further,  we  fired  a  shot  at  a  completely  different 
site  851  (see  Figure  1).  We  were  thus  able  to  evaluate  the  amplitude-yield  and 
distance  scaling  relations  found  from  the  KSG  series.  The  spectral  measurements 
at  Linac,  845,  and  858  gave  an  estimated  yield  of  94  ±  25  pounds  and  the  analog 
measurements  gave  50  t  20  pounds  for  a  100-pound  charge  at  Site  851.  This  rein¬ 
forces  the  conclusion  reached  in  the  preceding  paragraph. 

Amplitude  variations  associated  with  receiver  site  are  also  clearly  Iden¬ 
tifiable.  The  KSX  series  consisted  of  four  charges  detonated  at  802,  all  of 
nearly  equal  weight.  For  the  last  two,  KSX3  and  KSX4,  seismometers  at  Site  845 
were  moved  30  m  into  a  concrete  bunker  dug  Into  bedrock.  In  this  high-velocity 
material  we  expect  smaller  amplitudes  as  noted  by  other  investigators,  and  this 
effect  is  seen  in  Figure  5  for  the  KSX2  and  KSX3  shots.  The  respective  charge 
weights  of  these  shots  were  42  and  46  pounds  of  C4  equivalent.  The  ratios  of 
the  corresponding  "a",  "b",  and  spectral  average  measurements  were  1.49,  2.22, 
and  1.63,  showing  that  the  spectral  measurement  is  almost  as  good  as  "a"  and 
better  than  "b”  in  terms  of  rece 1 ver-s i te  independence. 


B.  Amplitude-Yield  Scaling  Relations 
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Table  2.  Ratio  of  Shot  X  to  KSF4 


X 

Seismic  spectral 

Seismic  analog 

expected 

comments 

KSF1 

.98  t. .35 

.52  4  .15 

1 

wet,  dark  mud  like  KSF4 

KSF2 

.48  i  .12 

.46  i  .07 

smal  1 

on  dry  sandstone 

KSF3 

.58  Z  .05 

.60  *  .05 

smal  1 

near  KSF2 

KSF5 

.85  1  .40 

.72  *  .04 

near  1 

near  edge  of  gully-mud 

KSF6 

.99  t  .46 

.62  .14 

1 

similar  to  KSF1 

KSF7 

.83  -  .13 

.80  4  .06 

1 

similar  to  KSF1 

KSF8 

1.27  4  .57 

.95  »  .09 

1 

similar  to  KSF1 
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To  determine  amplitude-yield  scaling  relations,  the  KSG  series  was  fired  at 
Site  802.  The  charges  used  were  24,  40,  60,  72,  93,  112,  167,  217,  316,  and  616 
pounds  of  C4  equivalent  high  explosive.  Regression  of  the  log  of  the  measure¬ 
ment  made  versus  the  log  of  the  charge  weight  in  pounds  was  performed  at  the 
three  recording  sites.  Results  appear  In  Table  3.  Values  of  the  yield  scaling 
exponent  “m"  range  from  0.6  to  1.0.  Tables  2  and  3  of  Murphy  and  Lahoud  (1969) 
give  slightly  higher  yield  scaling  exponents  for  displacement  than  for  velocity 
(0.761  versus  0.646  in  their  Table  3).  Here  the  spectral  measurements  represent 
displacement  and  the  analog  ones  represent  velocity.  Our  Table  3  gives  0.887 
and  0.774  for  the  respective  yield-scaling  exponents.  We  might  have  expected 
the  analog  measurements  to  scale  like  Springer  and  Hannon's  (1973)  Pn  data, 
which  gave  0.63  (similar  mode  of  wave  propagation).  As  found  by  those  authors, 
the  peak  frequency  in  the  analog  record  scales  weakly  with  yield,  here  varying 
by  less  than  25%  over  the  whole  range  of  charges  fired.  Undoubtedly  the  seismic 
waves  suffer  anelastic  attenuation  over  these  shallow  paths. 

DISCUSSION 

The  most  significant  result  found  is  that  spectral  measurements,  made  on 
wideband  digital  systems,  are  superior  to  analog  ones  In  several  ways  for  ampli¬ 
tude-yield  scaling.  First,  they  are  less  dependent  on  variations  in  source  cou¬ 
pling  and  receiver  conditions;  second,  they  permit  measurement 
standardization,;  third,  they  appear  more  objective.  We  discuss  each  of  these 
points  below. 


We  should  have  expected  less  dependence  on  scurce/rece  i  ver  effects  in  the 
spectral  averages.  The  analog  "a"  and  “b"  measurements  sample  only  the  first 
arrival,  a  particular  takeoff  angle  from  the  source.  The  spectral  measurements 
included  the  whole  group,  thus  samples  a  range  of  source  takeoff  angles.  This 
should  stabilize  the  measurement  somewhat  by  the  fact  that  it  is  an  average. 
Similarly,  the  receiver  response  function,  which  is  probably  strongly 
f requency-dependent ,  will  affect  the  measurement  less  if  we  average  over  a  size- 
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Table  3.  Yield-scaling  exponents,  log  A  =  m  log  Y  +  b,  where  "A"  is  amplitude, 
"Y"  is  yield  in  pounds. 


Measurement* 

m 

uncertainty** 

b  uncertainty** 

Linac  seismic  Z  6-9  Hz 

.937 

.0257 

-4.00 

.0524 

Linac  seismic  R  6-9  Hz 

.944 

.0339 

-3.81 

.0691 

Linac  seismic  T  6-9  Hz 

.798 

.0316 

-3.79 

.0645 

Linac  acoustic  Z  3-7  Hz 

.635 

.0539 

-2.54 

.110 

Linac  acoustic  R  3-7  Hz 

.710 

.0665 

-2.86 

.136 

Linac  acoustic  T  3-7  Hz 

.750 

.0585 

-2.61 

.119 

Linac  seismic  Z  "a" 

.782 

.0339 

1 .41 

.0691 

Linac  seismic  Z  "b" 

.725 

.0303 

1.98 

.0617 

Linac  seismic  R  "a" 

.809 

.0693 

1.06 

.147 

Linac  seismic  R  "b" 

.881 

.0360 

1 .46 

.0735 

858  seismic  Z  4-9  Hz 

.862 

.0363 

-4.00 

.0728 

858  seismic  R  4-9  Hz 

.880 

.0315 

-4.39 

.132 

858  seismic  T  4-9  Hz 

.798 

.0606 

-4.06 

.135 

858  acoustic  Z  1.5-3  Hz 

.735 

.0440 

-2.87 

.0883 

858  acoustic  R  1.5-3  Hz 

.925 

.0791 

-3.20 

.159 

858  acoustic  T  1.5-3  Hz 

.908 

.0646 

-3.36 

.130 

858  seismic  Z  "a" 

.815 

.0954 

1 .31 

.192 

858  seismic  Z  "b" 

.725 

.0575 

1 .89 

.115 

845  seismic  Z  8-12  Hz 

.695 

.0483 

-3.06 

.0985 

845  seismic  R  9-12  Hz 

.820 

.0661 

-2.76 

.135 

845  seismic  T  9-12  Hz 

1.25 

.189 

-3.73 

.386 

845  acoustic  Z  5-10  Hz 

.847 

.0738 

-2.84 

.150 

845  acoustic  R  5-10  Hz 

.659 

.0868 

-2.46 

.177 

845  seismic  Z  "a" 

.608 

.0441 

.2.46 

.0900 

845  seismic  Z  "b" 

.590 

.0438 

2.96 

.0892 

845  seismic  R  "a" 

.962 

.0578 

1  .65 

.118 

845  seismic  R  "b" 

.944 

.0449 

2.23 

.0915 

♦spectral  measurements  include  pass 

band  of 

average  used. 

In  this 

case,  the 

amplitude  is  in  units  of  micron-s. 

Analog 

"a"  and  "b"  measurements 

are  in 

arbitrary  units. 

**uncerta inties  of  slopes  and  intercepts  estimated  from  standard  regression 
formulas . 


able  pass  band  to  smooth  the  variations. 

Measurement  standardization  Is  made  possible  by  presenting  each  spectrum  as 
in  Figures  4  and  5,  with  ordinate  In  terms  of  absolute  ground  motion.  This  pei — 
mits  the  mixing  of  data  recorded  on  unmatched  systems.  It  also  eliminates  the 
problem  that  analog  measurements ,  taken  from  a  wideband  signal  viewed  through  a 
narrow-band  system,  cannot  be  related  simply  to  absolute  ground  motion;  thus, 
analog  measurements  from  different  systems  cannot  be  compared  directly  as  can 
absolute  spectral  ones. 

Objectivity  Is  aided  in  computing  spectral  averages  because  It  bypasses  the 
variation  in  signal  character  seen  on  analog  traces  by  discarding  the  phase  In¬ 
formation.  In  Figure  5  note  the  large  variation  in  signal  character,  and  In 
particular  the  large  difference  between  the  ratio  of  “a"  to  "b“;  however,  a 
spectral  average  basically  measures  the  power  In  the  whole  wave  group,  which 
should  be  a  more  stable  measure  of  yield  than  measuring  a  particular  phase  am¬ 
plitude,  whose  character  may  change  markedly  in  different  source  regions  as  Is 
seen  for  NTS  explosions. 

Finally,  note  that  the  wave  group  analyzed  was  a  refracting  phase  whose 
mode  of  propagation  is  evidently  similar  to  the  Pn  phase.  Thus,  our  results 
suggest  that  Springer  and  Hannon's  (1973)  amplitude-yield  scaling  relations  can 
profitably  be  extended  using  spectral  averages. 
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FIGURE  CAPTIONS 

Figure  1.  Left!  plan  view  of  standard  802  shotpoint,  special  851  shotpoint, 
and  the  three  seismic  recording  sites  described  herein.  Right!  plan  view 
of  shot  locations  in  the  KSF  source  coupling  experiment.  The  shotpoint 
numbered  ”4"  .is  the  standard  802  shotpoint. 

Figure  2.  Three-component  records  of  Shot  KSG7  as  recorded  at  the  Linac  site. 
Traces  from  top  to  bottom!  WWVB  time  code,  vertical,  radial,  and 
transverse  system  output,  essentially  ground  velocity.  All  seismic  traces 

are  scaled  to  the  same  plot  amplitude;  this  maximum  corresponds  to  the  am¬ 

plitude  given  left  of  each  seismic  trace  in  digital  counts. 

Figure  3,  Travel-time  curve  for  the  seismic  wave  as  recorded  for  all  shots  and 

all  distance  ranges.  A  least-squares  line  for  travel-time  "T"  in  seconds 

versus  distance  "X"  in  km  gives  T  =  0.327X  +  0.160. 

Figure  4.  After  windowing,  detrending,  and  tapering  (see  text),  the 
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vertical-component  seismic  phase  from  Figure  2  looks  like  the  lower  trace 
given  (maximum  amplitude  In  digital  counts  right  of  the  trace).  Above  In 
bold  Is  the  Fourier  amplitude  spectrum  after  smoothing  and  correction  to 
ground  motion.  The  thin  trace  Is  the  spectrum  of  record  segment  preceding 
event  onset  which  has  been  processed  the  same  way.  The  analog  measurements 
"a"  and  "b"  were  measured  as  Indicated  here. 

Figure  5.  Same  format  as  Figure  4,  but  as  recorded  for  KSX2  (left)  and  KSX3 
(right).  The  recording  site  for  KSX2  was  In  loose,  wet  soil,  and  for  KSX3 
it  was  30  m  SE  inside  a  concrete  bunker  set  In  bedrock. 


[53  THE  PN-RAYLEIGH  AND  ANOTHER  NEAR-REGIONAL  DISCRIMINANT 

WILLIAM  A.  PEPPIN 


SUMMARY.  A  neai — regional  seismic  discriminant  between  earthquakes  and  under¬ 
ground  explosions,  first  proposed  by  McEvIlly  and  Peppln  (1972)  and  based  on  a 
comparison  of  Pn  and  Rayleigh  amplitudes,  falls  for  a  selection  of  events  re¬ 
corded  at  a  new  wideband  station  Jamestown,  450  km  from  Nevada  Test  Site.  The 
same  events  can  be  correctly  identified  in  49  out  of  60  cases  If  a  comparison  is 
made  between  the  amplitude  ratios  of  envelopes  drawn  around  the  Sg  and  Pg 
phases.  This  discriminant,  similar  to  some  of  the  earliest  proposed,  has  the 
advantages  of:  potential  extension  to  small  (Richter  magnitude  ML  .It.  3.0) 
events  at  regional  distances;  ease  of  application;  ready  explanation  of  why  it 
is  effective;  and  probable  weak  dependence  on  earthquake  focal  mechanism. 
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McEvilly  and  Peppln  (1972)  and  Pepptn  and  McEvilly  (1974)  found  a  success¬ 
ful  seismic  discriminant  between  small  (ML  .ge.  3.6)  neai — regional  explosions 
and  earthquakes.  Using  wideband  (0.04  to  10  hz)  vertical  data  recorded  at 
Berkeley  (BRK),  Mina  (MNV) ,  Kanab  (KAN),  Landers  (LAN)  and  Elko  (ELK),  200  to 
550  km  from  Nevada  Test  Site  (NTS),  0.5  to  5-Hz  Pn  amplitude  was  plotted  versus 
0.02  to  0.1-Hz  Rayleigh  amplitude  for  explosions  and  earthquakes  (many  on  NTS) 
with  satisfactory  separation  of  explosions  from  other  events  at  all  stations  and 
all  magnitudes  studied  (3.5  .le.  ML  .le.  6.3).  In  1975  the  University  of  Cal¬ 
ifornia  at  Berkeley  began  operation  of  another  wideband  system  at  Jamestown 
(JAS)  in  the  Sierra  foothills.  An  attempt  was  made  to  check  the  effectiveness 
of  the  discriminant  based  on  new  (1975-1 979)  data  In  view  of  recent  Interest  In 
regional  seismic  discrimination.  I  present  results  of  the  analysis  for  this  and 
one  other  discriminant  in  this  paper. 

DATA 


In  Figure  1  is  presented  the  Instrument  characteristics  for  data  traces  re¬ 
corded  at  JAS.  Multigain  recording  of  the  4  data  channels  Is  on  compensated, 
40-percent  deviated,  FM  analog  magnetic  tape.  Data  channels  consist  of  a 
high-gain  short-period  Genioff  vertical  (SPZ),  a  wideband  velocity  system 
(0.025-10  Hz,  BBV),  and  a  wideband  displacement  system  (capacitance  transducer 
on  the  same  pendulum).  The  high-gain  and  low-gain  displacement  outputs,  DHG  and 
DLG,  are  separated  by  46  dB.  The  DLG  channel  Is  such  that  a  one  megaton  explo¬ 
sion  on  NTS  records  near  full-scale.  Therefore,  essentially  any  seismic  event 
in  California  or  Nevada  with  3.0  .le.  ML  .le.  6.5  Is  well-recorded  by  the  sys¬ 
tem  on  at  lease  one  data  trace,  providing  if  anything  superior  data  to  that  used 
in  the  earlier  studies.  Analog  playouts  of  filtered  Pn  (0.5  to  5  Hz,  48 
dB/octave)  and  Rayleigh  waves  (0.02  to  0.1  Hz,  48  dB/octave)  were  prepared  using 
an  experimental  procedure  identical  to  that  employed  by  McEvilly  and  Peppln 
(1972)  and  Peppin  and  McEvilly  (1974). 


SELECTION  OF  EVENTS  FOR  ANALYSIS 
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We  face  the  frustration  that  no  earthquake  has  occurred  on  NTS  between' 1975 
and  1979  that  would  be  suitable  for  direct  comparison  with  explosions  as  done  by 
Peppin  and  McEvIlly  (1974).  Thus,  following  McEvilly  and  Peppin  (1972),  we  take 
earthquakes  at  an  approximate  distance  to  JAS  equal  to  the  distance  to  NTS. 
Events  selected  are  presented  in  Tables  1  and  2. 

RESULTS:  PN  VERSUS  RAYLEIGH  DISCRIMINANT 

Figure  2  is  a  plot  of  Pn  versus  Rayleigh  amplitudes  for  the  events  given  In 
Table  1.  DHG  and  SPZ  were  taken  as  the  standards  for  measuring  Rayleigh  and  Pn 
amplitudes,  respectively.  Due  to  a  wide  range  of  amplitudes,  empirical  rela¬ 
tions  were  developed  between  these  standards  and  the  other  traces  available: 

DHG  =  BBV/<11.1  +  2.6) 
for  Rayleigh  waves  (22  common-  events) 

SPZ  =  BBVC17.1  ±  6) 
for  Pn  (15  common  events) 

SPZ  =  DHG(34.3  ±  12) 
for  Pn  (8  common  events) 

Using  these  and  the  static  gain  factor  of  200  between  DHG  and  DLG,  all  Rayleigh 
amplitudes  were  converted  to  DHG  and  all  Pn  amplitudes  were  converted  to  SPZ, 
the  numbers  given  in  Table  1.  Starred  readings  in  this  table  were  converted  to 
the  standard  ones  using  the  above  empirical  relations.  Figure  3  shows  the  loca¬ 
tions  of  JAS,  NTS,  and  the  earthquakes  used  in  this  study. 

RESULTS:  PG  VERSUS  SG  ENVELOPE 

The  appearance  of  the  records  under  considerations  suggested  another 
promising  discriminant.  In  Figure  4  an  envelope  was  drawn  around  the  Pg  and  Sg 
phases  on  visicorder  playouts.  The  amplitudes  A  and  B  were  read  and  the  ratio 
B/A  formed.  The  SPZ  channel  was  selected  as  the  standard  one  since  It  recorded 
events  of  smallest  magnitude.  Empirical  relations  were  developed  between  SPZ 
and  the  other  channels  as  follows: 

SPZ  =  DHG( 17.6  ±  7.8) 
for  P,  23  common  events 
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Table  1.  Amplitude  data  for  the  Pn-Rayleigh  discrimination,  measured 
in  arbitrary  units. 


IDENTIFIER 

EVENT  NO 

PN  AMP 

RAYLEIGH  AMP 

COMMENTS 

OW 

- 

975  JUN 

27  0727  J 

1  ’ 

32 

88* 

COLLAPSE 

975  JUN 

28  2948* 

2 

28 

77* 

COt  LAPSE 

975  JUL 

01 

0451 

3 

18 

30 

COLLAPSE 

975  JUL 

fll 

1814 

4 

40 

248* 

COLLAPSE 

975  SEP 

06 

1700 

5 

104 

18 

EXPLOSION 

i/b  OCT 

24 

1712 

6 

323* 

57 

EXPLOSION 

9/5  OCT 

28 

1430 

7 

170U0* 

16800 

EXPLOSION 

J/b  NOV 

20 

1500 

8 

2460 

3440 

EXPLOSION 

175  NOV 

26 

1530 

9 

104 

47 

EXPLOSION 

J/5  JAN 

03 

1915 

10 

6420* 

25600 

EXPLOSION 

976  JAN 

04 

0118 

11 

42 

76 

COLLAPSE 

9/5  JAN 

S4 

1616 

12 

18+ 

60+ 

COLLAPSE 

975  JAN 

17 

2140 

13 

22 

40 

COLLAPSE 

975  FEB 

04 

1420 

14 

12/5* 

1160* 

EXPLOSION 

9/6  FEIJ 

04 

1440 

15 

1500 

1309* 

EXPLOSION 

9/6  FEB 

12 

1445 

16 

5903 

5525 

EXPLOSION 

9/6  FEB 

14 

1130 

17 

3434* 

24000 

EXPLOSION 

976  FEB 

19 

1701 

18 

28 

10 

EARTHQUAKE.  AZ 

125,  MAG 

4.3.  D  332 

9/6  MAR 

DS 

1400 

19 

3485* 

6400 

EXPLOSION 

9/6  MAR 

14 

1430 

20 

7480* 

17600 

EXPLOSION 

9/6  MAR 

14 

1625 

21 

56 

200 

COLLAPSE 

9/5  MAY 

12.  1950 

22 

144 

88 

EXPLOSION 

9/6  JUN 

07 

0335 

23 

2 

44 

EARTHQUAKE ,  AZ 

115.  MAG 

4.1.  D  392 

378  AUG 

26 

1430 

24 

438 

248 

EXPLOSION 

9/6  DEC 

08 

1450 

25 

136 

112 

EXPLOSION 

•9/6  DEC 

21 

1539 

26 

104 

1 

EXPLOSION 

977  JAM 

13 

19/19 

27 

74 

10* 

EARTHQUAKE .  AZ 

165,  MAG 

3.8.  D  3/0 

97/  JAN 

13 

2009 

28 

44+ 

28 

EARTHQUAKE ,  AZ 

340.  MAG 

3.7.  D  3/7 

9/7  FEB 

16 

1/53 

29 

96 

22* 

EXPLOSION 

9/7  AUG 

04 

1640 

30 

328 

440 

EXPLOSION 

977  AUG 

12 

0220 

31 

224 

180 

EARTHQUAKE,  AZ 

165,  MAG 

4.8.  D  421 

977  AUG 

19 

1765 

32 

744 

1386* 

EXPLOSION 

I'/  SEP 

15 

1436 

33 

72 

11* 

EXPLOSION 

i//  ocr 

26 

1415 

34 

Cu 

76 

EXPLOSION 

•9/7  NOV 

10 

0235 

35 

48+ 

15 

EARTHQUAKE,  AZ 

145,  MAG 

4.0.  0  340 

9/7  NUV 

17 

1930 

36 

88 

26 

EXPLOSION 

)//  DEC 

14 

1530 

37 

3264* 

1400 

EXPLOSION 

9/8  FEB 

14 

0435 

38 

3'-*8 

256 

EARTHQUAKE.  AZ  075,  MAG 

4.8,  D  365 

;/U  MAY 

23 

0547 

39 

104 

74 

EARTHQUAKE ,  AZ  045,  MAG 

4.6.  0  434 

;/u  jut 

16 

5421 

40 

64 

15 

EARTHQUAKE ,  AZ 

166.  MAG 

4.3.  0  320 

.1/8  AUG 

01 

0902 

41 

56 

700 

EARTHQUAKE.  AZ 

346.  MAG 

4.6.  0  337 

AUG 

Ill 

034' 

42 

30 

528 

EARTHQUAKE.  AZ 

345.  MAG 

4.5.  f)  337 

:.'U  AUG 

111 

1025 

43 

44 

58 

EARTHQUAKE.  AZ 

34b.  MAS 

4.2.  0  337 

.'8  AUG 

31 

1400 

44 

5440* 

1400 

EXPLOSION 

9/5  AUG 

31 

2356 

45 

42 

86 

COLLAPSE 

!/:9  MOV  02 

1525 

46 

158 

65 

EXPLOSION 

'  JAN 

O'  20 

47 

3IJ+ 

39 

EARTH!  '.v\KE,  AZ 

75.  MA.G  ■ 

4.2,  0  365 

'  J.V 

,"l 

•!!!’  0 

48 

2.  •  J 

65  * 

tXI'LOSiO.K 

nrAQ'.NGS  BY  Co: .'VERSION  TO  STANDARD  ON!  S  (SEE  TEXT ) 

i*t  .3  p.e/.w.gs 
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SP2  =  DHG(13.6  ±6.1) 
for  S,  23  common  events 

SP2  -  BBV< 13.9  ±  4.2) 
for  P,  14  common  events 

SP2  =  BBVC13.0  ±  4.5) 
for  S,  24  common  events. 

In  Table  2  we  present  the  events  used,  the  P  amplitude,  and  the  S  to  P 
ratio  which  is  to  form  the  discriminant.  Starred  readings  In  the  Table  denote 
conversions  to  the  standard  ones  using  (2).  In  Figure  5  Is  plotted  S  to  P  ratio 
versus  log  of  P  amplitude. 

DISCUSSION:  PN  VERSUS  RAYLEIGH  DISCRIMINANT 

Figure  2  documents  the  complete  failure  of  the  Pn  versus  Rayleigh  discrimi¬ 
nant  for  the  JAS  data  set.  Note  In  Figure  3  that  BRK,  JAS,  and  MNV  are  all  lo¬ 
cated  NW  of  the  test  site.  Discrimination  based  on  Pn  and  Rayleigh  waves  was 
very  effective  at  BRK  and  MNV,  but  nonexistent  at  JAS  between  these  two  sta¬ 
tions.  As  Peppln  (1974)  noted:  the  theoretical  basis  for  such  a  discriminant 
hasn't  been  unequivocally  established,  so  perhaps  the  result  should  have  been 
expected.  » 

DISCUSSION:  PG  VERSUS  SG  ENVELOPE  DISCRIMINANT 

The  discriminant  shown  In  Figure  5,  sir.  i  ,ar  to  those  based  on  the  Lg  phase, 
appears  to  have  promise  and  could  certainly  be  refined  following  procedures 
adopted  by,  for  example,  Murphy  and:  Lahoud  (1975).  The  dashed  line  In  Figure  5 
is  a  decision  line  (explosions  below  and  earthquakes  above  with  one  exception, 
#28  in  Table  2).  The  dot-dashed  lines  are  based  on  estimated  uncertainties  of 
20  events  for  which  the  ratio  could  be  obtained  on  all  three  channels.  Thus, 
events  Inside  the  band  delimited  by  these  lines  are  not  identified  with  certain¬ 
ty.  Of  60  events,  11  are  unresolved  and  one  earthquake  is  misldentlf led  as  an 
explosion.  Collapses,  which  overlap  the  two  populations,  are  presumably  dis¬ 
tinctive  enough  as  to  cause  no  problem  for  Identification  by  other  methods.  No 


Table  2.  Data  for  the  discriminant  based  on  the  ratio  of  S  to  P  envelope 
amplitude. 


Identifier  Event  No.  P  Displacement  S/P  ratio  Comments 


975  OUN  27  0727 

1 

15. 

1.12* 

COLLAPSE 

375  JUN  28  0948 

2 

14. 

1.10* 

COLLAPSE 

975  OUL  81  D451 

3 

12. 

1.05 

COI.LAPSE 

1/5  OUL  01  1814 

4 

39. 

2.68 

COLLAPSE 

3/5  SEP  86  1/00 

6 

22. 

1.11 

EXPLOSION 

375  OCT  24  1712 

6 

67. 

1.24 

EXPLOSION 

3/5  OCT  28  1430 

7 

76. 

0.32 

EXPLOSION 

375  NOV  18  1530 

8 

12. 

1.75 

EXPLOSION 

375  NOV  20  1600 

9 

1400. 

0.50 

EXPLOSION 

3/5  NOV  26  1530 

10 

36. 

1.47 

EXPLOSION 

3/6  JAN  03  1915 

11 

4200. 

0.48 

EXPLOSION 

376  JAN  04  0118 

12 

13. 

1.26 

COLLAPSE 

376  JAN  04  1616 

13 

8. 

2.40 

COLLAPSE 

376  JAN  18  0/20 

14 

14. 

1.75 

COLLARS’7 

3/6  JAN  21  1840 

15 

9. 

2.11 

COLLAPSi. 

3/6  FEI3  03  0015 

16 

8. 

1.33 

COLLAPSE 

3/6  Ftl3  04  1420 

17 

480. 

1.20 

EXPLOSION 

376  FEB  04  1440 

18 

416. 

1.58 

EXPLOSION 

376  FEB  12  1445 

19 

1400. 

0.71 

EXPLOSION 

3/6  FEB  14  1130 

20 

3800. 

0.53 

EXPLOSION 

3/6  FEB  19  1701 

21 

16. 

3.97 

EARTHQUAKE,  AZ  125,  MAG 

4.3,  D 

332 

376  FEB  26  1450 

22 

8. 

1.28 

EXPLOSION 

376  MAP  09  1400 

23 

2000. 

1.33 

EXPLOSION 

176  MAR  09  1655 

24 

16. 

1.26 

COLLAPSE 

376  MAR  14  1450 

25 

4400. 

0.45 

EXPLOSION 

376  MAR  14  1525 

26 

30. 

1.22* 

COLLAPSE 

3/6  MAY  12  1950 

27 

63. 

1.18 

EXPLOSION 

376  JUN  87  0335 

28 

13. 

1.77 

EARTHQUAKE,  AZ  115,  MAG 

4.1,  0 

392 

3/6  JUN  19  1025 

29 

7. 

1.71 

EXPLOSION 

3/6  AUG  22  1014 

30 

6. 

2.2S 

EARTHQUAKE,  AZ  075,  MAG 

4.0,  D 

365 

3/6  AUG  26  1430 

31 

143. 

1.47 

EXPLOSION 

1/6  DEC  88  1450 

32 

69. 

1.12 

EXPLOSION 

376  DEC  21  1509 

33 

31. 

0.98 

EXPLOSION 

17/  JAN  13  0/19 

34 

13. 

2.00* 

EARTHQUAKE,  AZ  165,  MAG 

3.8.  D 

370 

1/7  JAN  13  2039 

36 

6. 

3.95 

EARTHQUAKE,  AZ  343,  MAG 

3.7,  D 

377 

i //  FEB  16  1753 

36 

25. 

1.25 

EXPLOSION 

37/  MAY  31  1640 

37 

3. 

3.67 

EARTHQUAKE,  AZ  330.  MAG 

3.7,  D 

330 

3/7  AUG  04  1640 

38 

130. 

0.86* 

EXPLOSION 

3/7  AUG  12  0220 

39 

40. 

3.73 

EARTHQUAKE.  AZ  165,  MAG 

4.8,  D 

421 

3/7  AUG  19  1/55 

40 

471. 

1.27* 

EXPLOSION 

3 77  SEP  15  1436 

41 

18. 

1.28 

EXPLOSION 

37/  OCT  26  1415 

42 

31. 

1.54 

EXPLOSION 

3 77  NOV  10  0236 

43 

14. 

3.28 

EARTHQUAKE.  AZ  145.  MAG 

4.0,  0 

340 

37/  NOV  17  1930 

44 

43. 

1.18 

EXPLOSION 

3 77  DEC  14  1530 

46 

49S. 

1.13 

EXPLOSION 

3/8  FEB  14  0435 

46 

15. 

2.96 

EARTHQUAKE.  AZ  075,  MAG 

4.8.  D 

365 

378  MAY  23  0547 

47 

20. 

4.45 

EARTHQUAKE,  AZ  045.  MAG 

4.6,  D 

434 

!'/8  JUl!  16  0421 

48 

64. 

2.17 

EARTHQUAKE,  A2  165,  MAG 

4.3,  D 

370 

378  J'JL  07  1480 

49 

24. 

1.16* 

EXPLOSION 

1/8  JUL  29  2232 

60 

62. 

2.08 

EARTHQUAKE,  AZ  058,  MAG 

4.0.  0 

380 

i/8  AUG  01  0902 

61 

40. 

3.87 

EARTHQUAKE.  AZ  345.  MAG 

4.6,  U 

397 

i/o  AUG  81  8947 

52 

24. 

5.27 

EARTHQUAKE,  AZ  345,  MAG 

4.5,  D 

397 

378  AUG  01  1026 

63 

16. 

2.21 

EARTHQUAKE,  AZ  345.  MAG 

4.2,  D 

397 

3/8  AUG  31  1400 

54 

1680. 

0.67 

EXPLOSION 

378  AUG  31  2356 

55 

21. 

0.50 

COLLAPSE 

3/8  NOV  02  1525 

56 

32. 

0.71 

EXPLOSION 

3/9  JAN  86  0120 

57 

16. 

2,55 

EARTHQUAKE,  AZ  75,  MAG 

4.2.  D  365 

3/9  JAN  24  1800 

58 

26. 

1.39 

EXPLOSION 

/•.;  Ft'i  22  0/16 

69 

8. 

4.16 

EARTHQUAKE ,  AZ  006.  MAG 

3.5.  0 

210 

79  1  Ell  22  1667 

60 

780. 

2.54 

EAR  1  •  IOUAKE ,  iV.  BSD,  MAG 

b.D.  D 

210 

BY  CONVERSION  TO  JAMESTOWN  SHUKT  PERIOD  (SFF  TEXT) 
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explosion  Is  mlsldentlf led  as  an  earthquake. 

The  following  can  be  cited  as  distinct  advantages  of  such  a  discriminant: 
(1)  it  Is  effective  to  magnitudes  as  small  as  can  be  seen  on  a  high-gain, 
short-period  Instrument,  l.e.,  less  than  ML  3.0  at  regional  distances;  (2>  It 
is  simple  to  apply;  (3)  it  Is  reasonably  easy  to  find  an  explanation  as  to  why 
it  Is  effective;  and  (4)  the  quantity  measured  is  a  ratio  that  probably  depends 
weakly  on  focal  mechanism  (many  takeoff  angles  at  the  source  contribute  to  Pg 
and  Sg),  hence  can  be  carried  over  directly  to  other  regions  without  need  of  a 
regional  calibration.  The  discriminant  should  work  because  envelope  amplitudes 
should  be  a  measure  of  S  energy  and  P  energy  leaving  the  source.  As  a  shear 
dislocation  releases  about  10  times  more  energy  in  S  than  in  P  (Haskell,  1964), 
it  is  expected  to  have  significantly  higher  S  to  P  ratio  than  an  explosion  in  a 
plot  like  Figure  5.  This  is  the  same  reason  given  by  Douglas  et  al.  (1971)  for 
effectiveness  of  the  Ms:mb  discriminant. 

CONCLUSIONS 

The  results  presented  here  Imply  that  the  Pn  versus  Rayleigh  discriminant 
between  explosions  and  earthquakes,  proposed  first  by  McEvIlly  and  Peppin 
(1972),  can  no  longer  be  considered  effective.  A  discriminant  based  on  relative 
excitation  of  the  Pg  and  Sg  phases  is  fairly  effective,  could  doubtless  be  re¬ 
fined,  and  holds  the  promise  of  extension  to  the  lowest  magnitudes  possible 
(less  than  ML  3.0  at  regional  distances).  However,  this  analysis  emphasizes  the 
need  for  theoretical  understanding  in  some  depth  before  empirical  seismic  dis- 
cr  i  m  i  .■ iants  can  be  used  confidently. 
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FIGURE  CAPTIONS 

Figure  1.  Instrument  characteristics  of  the  data  traces  recorded  on  analog 
magnetic  tape  for  the  Jamestown  wideband  system.  Curves  1  and  2 '•  DHG  and 
DLG;  3:  BBV;  4:  SPZ. 

Figure  2.  Log  of  Rayleigh  wave  amplitude  versus  log  of  Pn  amplitude.  The 
larger  symbols  denote  standard  measurements  (see  text)  and  the  smaller  ones 
correspond  to  any  of  the  starred  readings  from  Table  1.  Number  beside  each 
symbol  is  an  event  number  from  Table  1.  "Collapse"  events  are  the  hole 
collapses  following  many  of  the  explosions  on  NTS. 

Figure  3.  Map  view  of  the  Ca 1  If orn i a-Nevada  region  showing  locations  of 
wideband  station  JAS  and  NTS  (triangles),  and  the  earthquakes  used  In  this 
study  (numbers  correspond  to  table  2). 

Figure  4.  Sample  of  data  to  show  derivation  of  the  Sg-Pg  envelope  discrim¬ 
inant.  Envelopes  are  drawn  about  the  Pg  and  Sg  phases  by  hand,  and  the  am- 
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plltude  of  the  envelope  at  roughly  the  phase  onset  Is  measured,  the 
peak-to-peak  amplitudes  "A"  and  “B“  for  Pg  and  Sg,  respectively  In  the  fig¬ 
ure  . 

Figure  5.  Plot  of  Sg  to  Pg  amplitude  ratio  (Ratio  of  B  to  A  In  Figure  4> 
plotted  versus  log  of  Pg  envelope  amplitude,  the  values  given  In  Table  2. 
The  smaller  symbols  denote  starred  readings  from  the  Table.  Dashed  line  Is 
the  decision  line  (earthquakes  above  and  explosions  below).  Dot-dashed 
lines  delimit  a  region  of  uncertain  identification  (see  text). 


161  IMPULSE  SOLUTIONS  IN  A  LAYERED  ELASTIC  MEDIUM: 

EXTENDED  LAMB’S  PROBLEM 


Wi 11 iam  A.  Peppin 


ABSTRACT 

Johnson  (1974)  has  presented  the  complete  solution  for  displacements 
ar.d  scrasses  in  an  elastic  halfspace  resulting  from  a  localized  im¬ 
pulse  source  (Lamb,  1904).  This  paper  extends  Johnson’s  formulation 
to  a  multilayered  elastic  medium  by  presenting  a  solution  In  terms  of 
generalized  rays,  and  like  hi-  *h1s  method  is  exact  In  the  sense  that 
no  asymptotic  expansion  is  d.  Some  sample  calculations  are 
presented  which  compares  the  method  with  asymptotic  ray  theory  (Wig¬ 
gins  and  Helmberger,  1974),  and  which  shows  the  Importance  of  the 
near-field  terms  in  the  analysis  of  real  seismograms. 


1,  Introduction 


The  purpose  of  this  paper' Is  to  extend  Johnson's  (1974)  results 
for  the  elastic  displacements  resulting  from  point-force  and 
point-couple  sources  In  a  homogeneous,  Isotropic,  elastic  halfspace; 
the  methods  used  parallel  his  work  closely  (Sato,  1973a,  b  derives 
similar  results  using  another  approach).  The  motivation  for  this  pre¬ 
sentation  derives  from  recent  emergence  of  excellent,  wideband  digital 
data  close  to  explosion  and  earthquakes  sources  and  the  development  of 
the  seismic  moment-tensor  inversion  method  for  source  parameters  (Gil¬ 
bert,  1970;  Stump  and  Johnson,  1977;  Stump,  1979).  The  formulation 
presented  here  Is  not  amenable  to  fast  computation  of  such  solutions, 
and  the  reader  who  needs  speed  Is  referred  to  existing  generalized-ray 
codes  following  Wiggins  and  Helmberger  (1974)  or  to  the  very  fast  WKBJ 
code  written  by  C.  H.  Chapman,  J.  A.  Orcutt,  and  M.  Relchle  fol¬ 
lowing  Chapman  (1978)  and  Dey  Sarkar  and  Chapman  (1978).  It  finds  ap¬ 
plication  in  cases  where  impulse  solutions  are  required  very  close  to 
the  source  (within  10  km),  as  neai — field  effects  can  be  observed 
strongly  (see  Section  11).  Moreover,  in  many  cases  the  asymptotic  ray 
expansions  continue  to  be  accurate  even  at  close  distance  ranges. 
Thus,  a  second  motivation  is  to  explore  the  need  for  the  exact  formu¬ 
lation  at  these  ranges.  Only  aspects  of  the  mathematics  unique  to 
this  formulation  are  described  below;  the  interested  reader  will  find 
all  details  In  Peppln  (1974)  and  Johnson  (1974).  In  this  paper  equa¬ 
tion  number  "nnn"  from  the  latter  reference  Is  cited  as  "Jnnn". 

2.  Cnoice  of  Coordinates 

Let  a  right-handed  Cartesian  coordinate  system  (XI,  X2,  X3)  be 
imbedded  in  the  halfspace  with  the  plane  X3  =  0  coincident  with  the 
free  surface  and  with  the  X3-axis  passing  through  the  source,  which  is 
located  at  (0,  0,  X3'  =  -h).  Let  unprlmed  and  primed  coordinates  re¬ 
present,  respectively,  receiver  and  source  coordinates.  Note  that  X3 
is  positive  upward,  unlike  Johnson  (1974);  see  Figure  1. 


3.  Statement  of  the  Problem 
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Assume  a  halfspace  divided  Into  homogeneous,  Isotropic,  elastic 
layers.  These  are  numbered  downward  from  the  top  layer,  which  Is  1, 
and  Interfaces  are  numbered  downward  from  the  top  of  the  halfspace, 
which  is  a  free  surface  (interface  number  0).  Each  layer  is  defined 
by  the  four  quantities  a,  b,  r,  and  H,  the  P-  and  S-wave  velocities, 
the  density,  and  the  layer  thickness,  respect i vely.  All  numeric  sub¬ 
scripts  pertain  to  the  layer  number.  All  interfaces  are  assumed  weld¬ 
ed,  so  that  displacement  and  surface  traction  are  continuous  across 
them.  Assume  a  source  buried  at  depth  “h"  consisting  of  a  point  force 
in  the  direction  < f 1 ,  f2,  f3)  (cf  04): 

£ (*.t)  =  a.-fi,  -fs)T  — 

In  (1)  the  deltas  denote  Dirac  delta  functions,  the  source  is  located 
at  (0,  0,  -h),  and  the  source  time  history  is  No(t),  taken  zero  for 
all  times  t  <  0.  The  problem  is  to  compute  displacements  at  any  point 
on  the  surface  of  the  halfspace  in  response  to  this  source. 


4.  Notation 
h 

t — 

y  9 

( sec/cm  )2 
No(t ) 
p.q 

r 

£  & 
s 


source  depth  (cm) 

l  p— l  “t  r-> 

p2  +  q2  or  a2  +  )b,2,  squared  horizontal  wave  slowness 

time  dependence  of  applied  source 

wave  slownesses  in  x  and  y  directions  (sec/cm) 

epicentral  distance  (cm) 

Rayleigh's  function  (o2  +  412nanb)  -IT-  ■+■  A  A  t(  * 
temporal  Laplace  transform  variable  (sec)-l 


T 1 j 

i — ^ 

Cauchy  stress;  jth  component  of  traction  across  the 

plane  perpendicular  to  the  xi  axis  (dyne/cm2) 

ui 

elastic  displacement  (cm) 

oC 

P-wave  velocity  in  jth  layer  (cm/sec) 

P 

t— * 

S-wave  velocity  in  jth  layer  (cm/sec) 

(l/a^>2  -  ^J).5  vertical  P-wave  slowness  In  jth  layer 

( sec/cm)}  — 

— > 

X  t_j 
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7/5/ 

OK  <l)^j2  -  12). 5  vertical  S-wave  slowness  In  Jth  layer 

(sec/cm)  | - 9 

fjv^j  /Y  u  Lame's  elastic  constants  In  jth  layer  (dyne/cm) 

J  j  “j 

Vd  Pj  density  In  jth  layer  (g/cm3) 

(p  counterclockwise  azimuth  angle  from  E,  source  to 

receiver^ - '■> 

om^a  jCl  nb2  -  12  =  l/b2  -  212  rjf/  -  J/ p-  -til'1 

5.  Solutions  in  the  transform  domain  Following  Johnson  (1974)  we 
apply  coordinate  transformations  to  the  Cauchy-Navler  equations  J1  and 
the  isotropic  generalized  Hooke's  law  (Mase,  1970,  p  144)  assuming  In¬ 
finitesimal  strain  (Mase,  1970,  p  83).  These  coordinate  transforma¬ 
tions  are  defined  below;  they  take  time  “t" ,  “xl",  and  "x2“  into  "s", 
"p" ,  and  "q”,  respectively.  For  any  function  "f"  of  these  variables 
we  define 
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The  requirement  for  solution  of  the  problem  Is  that  all  of  the  trans¬ 
forms  defined  In  (2)  -  (4)  exist  and  that  the  inverse  transformations 
can  be  computed  on  some  contour  parallel  to  the  Imaginary  axis.  These 
conditions  are  met  under  08  (see  also  Peppln,  1974,  pp  157-158).  The 
result  Is  nine  equations  In  the  transformed  displacements  uj  and  trac- 

LJ 

tions  T3J .  These  are  rearranged  to  separate  P-SV  from  SH  terms  (Pep- 
pin,  1974,  pp  156-160  and  234-236).  The  result  was  presented  In  equi¬ 
valent  form  by  Gilbert  and  Backus  (1966)  and  Chapman  (1973): 
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In  (5)  VI  -  V4  are  the  transformed  displacement-stresses  associ¬ 
ated  with  the  P-SV  problem  and  V5,  V6  are  those  associated  with  the  SH 
problem  (see  Peppin,  1974,  237-239).  In  general,  this  separation  Is 
not  possible  in  the  (XI,  X2,  X3,  t)  domain. 

To  solve  (5)  for  V  the  method  is  to!  solve  above  and  below  the 
source;  set  conditions  joining  together  the  solutions  above  and  below 
the  source.  First,  note  that  solutions  of  (5>  can  be  taken  conveni¬ 
ently  as  a  linear  combination  of  the  eigenvectors  of  the  matrix  M. 
Chapman  (1973)  has  presented  these  and  the  corresponding  eigenvalues 
in  terms  of  vertical  wave  slownesses 


where  "c"  is  a  or  b. 
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where 


n  --  nr>  -A 

A  ~  nb2  -~12Q- 


/S'  A 


Peppin  (1974),  pp  240-245  shows  that  A,  B,  and  C  are  upward-travelling 
P,  SV,  and  SH  waves  while  D,  E,  and  F  are  downward-travelling  P,  SV, 
and  SH  waves.  Thus,  we  solve  (5)  above  the  source  In  terms  of  A,  B, 

.  'V'  .  x* 

and  C,  and  below  the  source  in  terms  of  D,  E,  and  F.  That  Is,  we  seek 

Ay  Ay 

constants  K1 ,  K2,  and  K3  such  that  above  the  source 
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and  constants  K4,  K5,  and  K6  such  that  below  the  source 

V  :K4oe  f  \c^Ee  -t-i4bhe  —  ( v; 


Substitution  of  (8)  and  (91  in  (5)  shows  that  these  are  valid  solu¬ 
tions  for  VI,  V3,  and  V5  for  any  value  of  the  constants  "K“.  To  find 

i— i  *-» 

these  constants,  note  that  V2,  V4,  and  V6,  the  transformed  stresses, 
suffer  discontinuities  due  to  the  presence  of  the  delta  function. 
Thus,  from  (5) 
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Thus,  inserting  (10)  in  (5)  to  match  the  solutions  above  and  below  the 
source 
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where 

H  ~ 

-  (  f ,  f  + 

f  2  S'  ) 

T  5 

-  (  f  2  S>  *  — 

f<  ^  s) 

z  = 

-  fj 

The  solutions  (8)  and  (9)  to  (5)  are  obtained  be  Inserting  (12) 
and  (13).  If  the  source  Is  on  an  Interface,  (12)  and  (13)  do  not 
apply,  and  we  would  have  to  return  to  (11)  and  solve  with  differing 
source  properties  above  and  below  the  source. 

6.  Solutions  at  the  Free  Surface 

All  solutions  presented  in  this  paper  assume  a  burled  source  and 
a  receiver  placed  on  the  free  surface.  The  condition  to  be  set  Is 

— r  w-tj  t: 

that  the  tractions  V2,  V4,  and  V6  vanish  there.  The  three  cases  to  ' 

«—  w-i  t_) 

consider  are:  Incident  P,  SV,  or  SH.  Notation  used  Is  given  In  Fig¬ 
ure  2.  With  these  definitions,  we  see  for  the  P  problem  that 

V  (  -  °)  -  £  ■**  D  r  E  - ( 

. — ■  ^ 

With  the  transformed  surface  tractions  V2  =  V4  =  0,  we  find  that 

p.  ■  7  _  ~  JL  -  -4  £  I] I'i  *  pq  •  ii  {  ftf‘7  x  <7;;)  z 

-  -  -  -•  - .  ° -t l  - 

*  ---nO 

(&  --  Jl2-  f 

For  the  incident  SV  problem  we  find  P022  =  P044,  P042  =  P024.  Thus, 

i — .  — •  _j  i — j 

from  (14)  and  similar  relations  for  SV  and  SH  we  find  that  the  Inci¬ 
dent  displacements  are  modified  by  Interaction  with  the  free  surface 
through  multiplicative  factors.  The  transformed  solutions  at  the  free 
surface  apart  from  the  exponential  factors  in  (8)  and  (9)  are 
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In  (16)  *Kj",  j  =  1  to  6  are  the  source  factors  from  (12)  and  (13); 

—i 

all  elastic  properties  and  vertical  slownesses  depend  upon  the  propel — 
ties  of  the  source  layer  in  these  expressions.  In  contrast,  all  elas¬ 
tic  properties  and  vertical  slownesses  appearing  in  (16)  apart  from 
the  Kjs  depend  on  the  properties  of  the  uppermost  layer.  Note  how 

i-J 

(16)  shows  complete  separation  of  the  effects  of  source  and  receiver 
in  this  problem. 

7.  Shift  to  Polar  Coordinates  in  the  Transform  Domain 


A  final  step  is  necessary  before  Inversion  to  the  desired  dis¬ 
placements.  Note  that  VI  and  V5  are  not  the  transforms  of  physically 

t-i  — * 

useful  quantities.*  see  (5).  To  recttfy  this,  DeHoop  (1961)  suggests 
a  shift  in  the  transform  parameters  as  follows: 
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where  "pX*"  Is  the  polar  angle  of  the  line  from  source  to  receiver 
measured  from  due  east  (Figure  1).  This  simplifies  the  Inversion  in¬ 
tegrals,  because  under  (17)  ,  > 

\/  f9 

-spxl  -sqx2  =  -sprcos(phi )  “  sqrsin(p)^l)  =  -sra 

and 

dpdq  =  dadb. 


Mow  define  longitudinal  and  transverse  components  in  terms  of  previ¬ 
ously-defined  quantities: 
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where  radially  away  from  the  source  and  transverse  to  the  right  look- 


ing  at  the  source  from 

the  receiver  are  taken  as  positive.  From  (17), 

(18),  and  (5) 
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At  this  point,  apart  from  the  modifications  needed  for  computing 
an  arbitrary  generalized  ray,  the  method  follows  Johnson  (1974)  at  and 
beyond  020.  The  only  significant  difference  up  to  this  point  Is  that 
the  SV  and  SH  waves  are  here  separated,  while  In  Johnson  (1974)  they 
are  not.  The  separation  is  required  in  the  multilayer  case  because  SV 
and  SH  ref lcct i on-transmi ss i on  coefficients  differ. 

7.  Inversion  to  the  Time  Domain 

We  present  here  specific  formulas  for  the  halfspace  case;  the 

modification  of  these  to  arbitrary  generalized  rays  Is  easy  (see 

below).  Applying  (18)  and  (19)  in  (8)  or  (9)  using  (12)  or  (12)  al- 

u«5 

lows  solution  in  the  transform  domain  for  (ur,  utl^ta,  uz).  To  these 

■  .  »-  j 

apply  inversion  integrals  (3)  and  (4).  Shift  the  Integration  vari¬ 
ables  from  dpdq  to  dadb.  The  result  at  the  free  surface  got  using 
(16)  is 
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.  A  /\ 

where  f  =  (fr,  fth^ta,  fz),  P  Is  similar  to  M  In  J16,  and  where  SV  + 

SH  is  similar  to  N  in  016,  i.e., 
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For  any  specified  source  type,  the  Integrals  (20)  must  be  reduced 

to  a  standard  form.  Here  we  show  the  reductions  expllctly  for  a 

source  of  type  (1).  The  integral  over  "a"  is  reduced  by  the  principle 

oc 

of  reflection  (Peppin,  1974,  pp  246-247)  to  the  domain  [0,iiVf).  Then 

&  ' 

we  shift  the  "b"  Integral  to  the  real  axis  (b  =  ip),  on  the  domain  [0, 
i)\f).  In  this  new  domain  we  discard  all  terms  in  the  integrand  that 


are  odd  in  X-  The  result  for  the  halfspace  case  is 
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where  P,  SV,  and  SH  are  derived  from  (20)  by  substituting  b  =  1e  and 
dropping  all  terms  odd  In  this  new  Integration  variable. 


Now  shift  the  Inner  "a"  Integrals  In  (22)  over  to  the  real  vari¬ 
able  "t"  which  has  units  of  time* 

t  £  c\r  +•  r/ch  t  c  -  *  or  p  -  —  m) 


This  deforms  the  contour  of  integration  from  the  Imaginary  axis  onto 
the  Cagnlard-deHoop  contour,  which  runs  along  the  real  "a"  axis,  then 

steps  into  the  first  quadrant  before  It  reaches  the  branch  point  at 

_ 

1  (etac)  =  0  on  the  real  “a"  axis.  The  deformation  requires  that  no  sin¬ 
gularity  lies  between  the  two  contours,  which  is  true  as  long  as  the 
source  is  below  the  free  surface  (Peppin,  1974,  pp  255-258).  Under 
this  variable  change,  the  first  integral  in  (22)  becomes 
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It  Is  crucial  for  the  method  to  note  that  in  (24)  the  lower  limit 
of  the  Inner  Integral  "tmln"  always  exceeds  zero.  This  Is  true  be¬ 
cause  for  a  range  of  values  t  >  0,  real  values  of  BaB  are  required  In 
(23),  so  that  the  Integrand  in  (24)  has  no  Imaginary  part. 

Physically,  Btmin"  corresponds  either  tc  the  Fermat  arrival  time  O 

(where  the  Cagniard-deHoop  contour  leaves  the  real  axis)  or  to  the 

earliest  headwave  arrival  time  (contour  encounters  a  branch  point 

along  the  real-"a"  axis  corresponding  to  a  zero  of  an  t)c  not  appearing  X 

explicitly  in  (23)).  For  the  halfspace  and  other  special  cases, 

"tmin"  can  be  given  explicitly  (see  J28).  Thus,  because  ’ t m 1 n"  is  po-  X 

sitive  in  all  cases,  when  in  (24)  we  invert  the  order  of  integration, 
the  inner  integral  (now  over  e)  extends  to  a  bounded  upper  limit.  A 
second  crucial  requirement  is  that  "tmln"  is  monotonically  increasing  * 

tr  - * 

in  "e",  or  else  this  order  reversal  cannot  be  accomplished  simply.  A  A 

"V 
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This  point  Is  explicitly  true  for  the  halfspace  case  <see  J27)  and  can 
be  proved  for  the  general  case  as  well  (see  below).  Under  this  furth- 


/-C i(t) 
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er  change,  (24)  becomes 


If  No(s>  Is  1/s,  then  No(t)  Is  H(t),  the  Heaviside  step  function.  In 

LJ 

this  case  the  factor  "s"  is  contained  only  in  the  exponent  of  “e*  In 
(25),  so  that  the  time-domain  results  we  want  consist  of  the  Integrand 
of  the  ”t“  integration  in  (25).  Thus,  the  time-domain  Inversion  of 
(20)  for  a  step-function  source  Is 
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The  form  (26)  Implies  that  the  displacements  can  be  computed  at 

any  specified  time  through  real  quadratures.  Once  the  desired  time 
"t”  Is  specified,  we  compute  the  upper  limit  e(t)  In  (26)s  see  below. 

fc 

Then  to  evaluate  the  Integrand  at  a  submesh  point  ^ " ,  we  solve  for 
"a”  in  (23).  All  expressions  within  the  integrand  are  functions  of 

fc 

"t",  "el",  and  the  value  of  "a1  just  found. 


8.  General Izat ion  to  Higher-Order  Sources 

The  displacements  resulting  from  highei — order  force  moments, 
equivalent  to  spatial  derivatives  of  the  above  fundamental  source  with 
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respect  to  source  coordinates,  can  be  obtained  as  suggested  by  Johnson 
(1974).  Suppose  that  we  differentiate  the  equations  to  be  transformed 
with  respect  to  source  coordinates  XI'  or  X2‘ .  Then  this  source  deri¬ 
vative  would  be  applied  to  (1).  Notice  that  an  Increase  In  source 
coordinates  Is  equivalent  to  a  decrease  In  receiver  coordinates  by  the 
same  amount.  Thus, 
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Under  (3)  and  (4)  the  unprtmed  (receiver)  XI  and  X2  derivatives  lead 

La  l~J 

to  a  multiplicative  factor  -sp  and  -sq,  respectively,  which  shows  up 
in  the  source  term  of  (5).  Thus,  the  primed  XI  and  X2  (source)  derl- 
vatives  multiplicative  factors  +sp  and  +sq  In  the  source  term  of  (5), 
so  that  derivatives  with  respect  to  either  source  or  receiver  coordi¬ 
nates  XI  and  X2  modify  the  solution  by  introducing  the  above  factors 

— *  U 

in  (21).  For  each  such  derivative  another  factor  of  ”s"  is 
introduced;  thus,  to  effect  Cagn i ard-deHoop  inversion  leading  to 
(26),  we  must  modify  No(t)  to  eliminate  any  accumulated  multiplicative 
factors  of  "s"  outside  the  exponent  in  (25). 

Consider  derivatives  with  respect  to  source  depth.  It  Is  conven¬ 
ient  to  carry  that  differentiation  through  to  (8)  and  (9)  and  take 
this  derivative  explicitly.  For  upgolng  rays,  noting  that  the  deriva¬ 
tive  vs  want  is 

d/d(-h), 

we  multiply  the  solution,  hence  (22)  by  +^jc  and  -s)*p  for  upgolng  and 

downgoing  rays,  respectively  (c  =  or.lj). 

(A  p 

Thus,  to  proceed  to  any  higher-order  source,  Introduce  In  (22) 
one  such  factor  for  each  derivative,  drop  all  terms  odd  in  "b"  from 
the  result,  and  proceed  exactly  as  above  to  a  real-valued  quadrature. 
Rather  than  present  a  list  of  these  solutions,  the  reader  can  see  most 
easily  how  these  derivatives  are  Introduced  by  studying  the  relevant 
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subroutine  for  computing  source,  receiver,  and  source-derl vatlve  fac¬ 
tors?  see  the  Appendix. 

9.  Modification  for  an  Arbitrary  Generalized  Raypath 

Formula  <26)  can  be  used  with  slight  modification  for  the  case  of 
an  arbitrary  raypath,  as  shown  In  Figure  3.  In  this  case  we  start 
with  the  Integral  for  SV  leaving  the  source  downward  (the  second  in¬ 
tegral  of  (26);  the  factor  K5^in  (12)).  Cisternas  et  al.  (1973)  j 1 

show  that  propagation  along  any  realizable  raypath  can  be  accounted 
for  in  the  transform  domain  by  products  involving  the 
ref lect i on-transmi ss i on  coefficients  for  each  Interface  interaction 
and  exponential  decay  factors  dependent  on  the  vertical  layer  thick¬ 
ness  traversed  (this  is  the  distance-decay  factor).  The  results  for 
the  ray  shown  In  Figure  3  is  to  replace  the  factor 

4-  fhj 

£ 

in  (9)  with 

-  >  '//s,  Ou  *-  -h)  -  S  'lh  Hz  ~  srl*,  H< 

$  (24  C  —  (Z7) 


where  S212f  and  5141  are  the  corresponding  reflection-transmission 

coefficients,  and  where  single  subscripts  pertain  to  layer  number. 

Thus  the  integrand  in  (26)  will  contain  the  added  multiplicative  fac- 
S  =  .i 

tor  (S2T2^t24)  and  (23)  explicitly  becomes 

j-  -  a  r  +  I'JSj.t  Hi  ~h)  '/  **■  f-  Hi  (2$) 


1 _ i  ^Zi-2. 

4  > 


Note  that  (23)  can  be  solved  analytically  for  "a"  (see  031  and 
032),  but  (28)  must  be  solved  numerically  for  "a".  For  real  values  of 
"a",  (28)  sometimes  has  two  real  roots;  the  smaller  of  these  is 
taken . 
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Thus,  the  procedure  is  to  specify  all  of  the  generalized  raypaths 
of  Interest,  introduce  the  corresponding  reflection-transmission  fac¬ 
tors  of  each  generalized  ray  In  the  integrand  of  (26)  and  proceed  to  a 
real-valued  quadrature  as  before  for  the  halfspace  case. 

Note  that  most  of  the  time  spent  In  numerical  evaluation  of  such 
integrals  is  in  finding  the  complex  roots  of  (28).  Also  note  that, 
for  complicated  raypaths,  the  kinematic  and  dynamic  groups  must  be  re¬ 
cognized,  as  a  single  member  of  such  a  group  has  no  physical  meaning. 

10.  Numerical  Aspects  of  the  Quadrature  in  (26) 


In  (26)  it  remains  to  discuss  the  computation  of  the  upper  limits 
of  integration.  These  must  be  computed  at  each  time  step.  The  topol¬ 
ogy  of  the  domain  of  integration  involving  these  limits  is  shown  sche¬ 
matically  in  Figure  4.  The  shaded  regions  in  the  figure  denote  re¬ 
gions  of  e~t  space  that  contribute  to  the  integral;  for  a  specified 
time  the  domain  of  integration  is  a  vertical  line  from  the  e  =  0  axis 
at  that  time  to  the  upper  limit  of  the  shaded  region.  Suppose  we  de¬ 
sire  to  compute  the  Integral  at  time  "tp"  (see  the  figure).  Then  the 
three  limits  el,  e2,  and  e3  must  all  be  known.  This  is  because  the 
integrand  in  (26)  contains:  bounded  but  infinite-slope  singularities 
at  el  and  e2;  a  square-root  singularity  at  e3.  To  achieve  sufficient 
accuracy  (1  part  in  500,000),  the  integral  is  subdivided  into  regions 
between  these  limits.  A  simple  change  of  variable  is  sufficient  to 
eliminate  each  part  of  the  Integrand  that  Is  non  polynomlc  (Peppln, 
1974,  pp  192-194). 


e,  u 


i_i  e 


The  headwave  limits  like  el  and  e2  derive  from  analytical  formu-  «-J  G 

l~l  L— » 

las  like  the  second  of  027  with  the  obvious  extension  to  the  multlla- 

JLU  1  /c*L 

yer  case  (i.e.,  solve  for  e  in  (28)  for  which  and  this  the  /  ' 

"  Am 

correspond i ng  headwave  slowness  is  less  than  any  of  the  slownesses  ap- 

•k;j-  ,  > 

pearing  in  the  nc^  of  (28)).  The  Fermat  limit  e3  can  only  be  given  uJ  ,A  ' 
explicitly  (first  of  027)  for  rays  in  which  only  one  vertical  slowness 
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nc  appears  In  (24),  such  as  internal  layer  reflections  In  addition  to 
the  halfspace  rays.  Otherwise  it  is  found  as  the  value  of  e  which  sa¬ 
tisfies  the  expression  corresponding  to  (28)  for  the  given  time  tp  and 
which  satisfies  the  saddlepoint  condition 


d  t 

cJ 


=  O 


Peppin  (1977)  shows  that  these  two  conditions  are  met  for  the  value  of 
“X“  satisfying 


~BU)  v  tt-  2  Hi'U  ){£  )- 

j  -  *  ;  -i  J  / 


-  o 


x  J 


where  the  sums  extend  over  each  distinct  vertical  slowness  in  expres¬ 
sions  like  (28)  for  a  particular  generalized  ray.  The  zero  of  fX^)  is 
easy  to  find  numerically;  only  one  real  root  exists  (Peppin,  1977). 
Substituting  this  root  In  (29)  permits  solution  for  ”a",  then  for  the 
desired  upper  limit  as  (^  -  12)1/2. 

a~  -j2* 


x  X 


n 


A  > 


c- 


A  second  aspect  of  the  numerical  calculations  Is  that  the  gener¬ 
alized  reflection-transmission  coefficients  must  be  known.  These  are 
provided  by  Chapman  (1973).  The  computation  of  these  can  be  labori¬ 
ous.  consuming  up  to  25%  of  the  total  processing  time  required.  To 
make  a  run  Involving  25  generalized  rays  in  a  fairly  complicated  layer 
stack,  100  seconds  of  CYBER176  time  are  required  for  close-in  calcula¬ 
tions  (at  distance  this  can  be  reduced  by  an  order  of  magnitude). 
Many  possibilities  exist  to  shorten  this  processing  time  considerably, 
but  these  have  not  been  fully  explored  because  of  the  existence  of 
much  faster,  approximate  methods  such  as  the  WKBO  method  of  Chapman 
(1978). 


11.  Numerical  Examples 

A  computer  program  MEXEC,  written  in  FORTRAN,  evaluates  (26)  for 
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any  generalized  raypath.  The  source  Is  specified  In  terms  of  an  arbi- 
trarl ly-or iented  point  force  <1)  or  of  spatial  derivatives  of  this 
source  with  respect  to  source  coordinates.  The  subroutine  appearing 
in  the  Appendix  Is  written  to  accommodate  extension  to  highei — order 
sources.  In  the  current  version  of  MEXEC  ten  components  of  ground  mo¬ 
tion  are  computed;  they  are  the  canonical  source  components  of  Stump 
and  Johnson  (1977)  and  Stump  (1979).  Any  point  source  can  be  ex¬ 
pressed  as  a  linear  combination  of  the  canonical  sources  through  the 
seismic  moment  tensor  (Stump,  1979);  these  sources  are  the  starting 
point  for  Stump  and  Johnson's  procedure  to  Invert  seismograms  for 
fault  parameters.  Output  of  the  code  has  been  checked  against  the 
calculations  presented  by  Johnson  (1974)  for  the  halfspace  and  by  a 
modification  of  0.  V.  Helmberger's  code  for  asymptotic  ray  calcula¬ 
tions  of  the  canonical  sources  (Walt  Silva,  personal  communication, 
1980).  MEXEC  runs  successfully  on  five  CYBER-based  systems  and  on  one 
32-bit  system. 

A.  The  Importance  of  Near-field  Motions 

In  Figure  5  are  presented  displacement  seismograms  recorded  on  a 
wideband  digital  system  (Peppin  and  Bufe,  1980)  close  to  the  hypo- 
centers  of  two  sizeable  earthquakes.  Also  shown  are  predicted  ground 
motions  computed  by  MEXEC  for  a  halfspace.  Note  In  the  figure  the 
distinct  ramplike  motion  between  the  P  and  S  phases.  This  Is  the 
near-field  motion  which  disappears  at  distance  from  the  source. 
Because  theoretical  and  observed  traces  are  so  similar,  little  doubt 
-emains  that  we  have  correctly  identified  near-field  motions  in  these 
cases.  Such  motion  has  been  observed  on  hundreds  of  displacement  se¬ 
ismograms  take  within  20  km  of  the  source  (Peppin  and  Somerville, 
1980). 


B.  Comparison  with  Asymptotic  Ray  Theory 


How  significant  an  error  is  made  by  the  utilization  of  the  much 
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faster  (factor  of  50)  asymptotic  ray  calculation?  In  many  cases,  al¬ 
most  none  as  we  shall  see.  A  test  run  was  made  for  a  source  depth  of 
8  km,  a  receiver  distance  of  32  km,  and  a  layer  thickness  of  32  km. 
The  results  agree  quite  closely  for  all  of  the  ten  canonical  source 
components  except  the  transverse  ones.  Figure  6  shows  a  comparison 
for  strlkeslip  sources  and  Figure  7  for  dipslip  ones.  The  disagree¬ 
ment  between  the  two  formulations  Is  significant  only  in  the 
near-field  part  of  the  transverse  components*  asymptotic  ray  theory 
produces  no  near-field  term.  The  data  in  Figure  5  require  these  terms 
to  explain  the  observed  ramp  between  the  P  and  S  phases.  Note, 
however,  that  the  step  offset  on  the  transverse  components  (the 
"far-f ield"  term)  agrees  quite  closely  in  spite  of  the  discrepancy 
suggested  by  the  figure  (for  drawing  convenience  only).  Also  note 
that  the  differences  between  the  two  theories  are  emphasized  in  that 
the  calculations  shown  are  for  a  ramp  function  source. 

C.  Source  Below  a  Layer 

Finally,  consider  the  case  of  a  source  below  a  slow,  soft  layer. 
In  Figure  8  we  compare  radial  and.  vertical  components  for  strlkeslip 
and  dipslip  sources.  Notorious  Is  that*  (1)  the  strlkeslip  source 
first  motion  Is  impulsive  and  the  dipslip  is  more  similar  to  the  he- 
adwave  onset  SPS,  seen  clearly  on  the  strlkeslip  components;  (2) 
first  arrivals  are  more  dominant  on  the  vertical  components  than  on 
the  radial  ones;  and  (3)  a  long-period  arrival,  probably  an  interface 
wave,  begins  about  4  1/2  seconds  after  the  origin  time.  Compare  the 
complexity  of  these  traces  with  the  explosion  source  of  Figure  9. 
Generally,  S  leaving  the  source  renders  the  seismograms  much  more  com¬ 
plicated.  The  calculations  for  the  ten  source  components  required  15 
seconds  of  CYBER  176  time  to  compute. 
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FIGURE  CAPTIONS 

Figure  1.  Illustrates  geometry  of  problem.  Source  is  at  "S",  at  a 
depth  "h"  below  the  free  surface  X3  =  0.  Receiver  is  at  epicentral 
distance  "r”.  ”phi“  is  measured  counterclockwise  from  east. 

Figure  2.  Gives  notation  for  three  problems  P,  SV,  or  SH  Incident  at 
free  surface.  The  factors  POij  are  the  reflection  coefficients  at  the 
free  surface.  Wavy  lines  denote  travelling  S-waves,  straight  lines 
denote  P-waves.  Vector  notation  is  defined  in  (7). 

Figure  3.  The  particular  ray  described  In  text.  Leaves  source  as  SV 
in  the  second  layer,  converts  to  P  on  reflection,  and  transmits  to 
free  surface  as  P.  Epicentral  distance  is  ”r",  source  depth  is  "h”. 

Figure  4.  Description  of  domain  of  integration  (shaded)  for  a  genet — 
all  zed  ray  such  as  shown  in  Figure  3,  with  two  headwaves,  arrival 
times  tHl  and  tH2,  and  a  Fermat  time  of  tF.  The  true  diagram  appears 
as  the  small  one  above  (each  front  monotonically  Increases  with  t). 
For  clarity  the  vertical  scale  is  reduced  to  an  appropriate  slowness 
and  expanded  in  the  larger  diagram  below.  Let  the  time  "tP"  be  se¬ 
lected  at  which  to  perform  the  quadrature.  Then  the  domain  of  integ¬ 
ration  is  the  vertical  strip  at  tP  covering  the  shaded  region.  As 
descr  ed  in  text,  we  subdivide  the  domain  into  three  separate  Inter¬ 
vals  using  the  limits  el,  e2,  and  e3.  Times  tC2  and  tCl  are  the  he- 
adwave  cutoff  times  beyond  which  the  headwave  fronts  have  no  physical 
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meaning.  The  Fermat  Integral  extends  to  all  times  (area  under  curve 
labelled  *tF"  In  top,  smaller  diagram). 

Figure  5.  Theoretical  versus  observed  displacements  for  two  earth¬ 
quakes  (left:  26  May  1980  1925  GCT,  Mo  2E22 . dyne-cm;  right:  26  May 
80  1857  GCT,  Mo  1.7E23  dyne-cm).  The  eplcentral  distance  and  the 
focal  depth  are  each  about  5  km  for  these  events.  For  the  smaller 
event  a  rise  time  of  0.15  second  Is  assumed,  and  for  the  larger  event 
a  rise  time  of  0.45  second.  The  observations  are  clipped  after  onset 
of  S  on  the  larger  event.  The  "neai — field  terms”  give  the  prominent 
ramp  between  the  P  and  S.  Synthetics  assume  a  halfspace. 

Figure  6.  Radial,  transverse,  and  vertical-component  strikesllp  seis¬ 
mograms  for  a  source  in  a  layer  32  km  thick.  Solid  traces:  this  for¬ 
mulation;  dotted  traces:  asymptotic  Cagniard  theory  (Wiggins  and 
Helmberger,  1974).  Rays  include  the  directs  (P,  S),  the  free-surface 
headwave  (sP),  and  the  first  multiples  (PP,  PS,  SP,  SS).  The  vertical 
bars  give  amplitudes  of  2E~25,  5E-25,  and  IE-25  cm,  respectively,  for 
a  source  with  a  strength  of  1  dyne-cm.  Source  time  function  is  a 
ramp. 


Figure  7.  Same  format  as  Figure  6  for  a  dlpsllp  source;  radial, 
transverse,  and  vertical  components.  The  vertical  bars  are  IE-25, 
IE-25,  and  2E-25  cm  for  a  source  as  in  Figure  6. 

Figure  8.  For  a  source  below  a  soft,  slow  layer  1  km  thick.  Left: 
strikeslip  radial  and  vertical;  right:  dlpsllp  radial  and  vertical. 
Raysum  Includes  all  directs  and  first  multiples  arriving  before  the  S 
wave.  "SPS"  is  the  headwave  associated  with  the  direct  ray  SS,  the 
S-wave.  The  source  depth  and  receiver  distance  are  1.3  and  8  km,  res¬ 
pectively.  the  layer  velo-ities  and  density  are  1.7E5,  1E5,  and  2.0 
cgs,  and  the  halfspace  velocities  and  density  are  3.34E5,  1.95E5,  and 
2.7  cgs.  On  each  trace  the  bar  represents  IE-23  cm  of  displacement 
for  a  source  with  strength  1  dyne-cm.  Source  time  function  is  a  step 
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function,  and  tha  Haight  of  tha  splkas  ara  conditioned  by  the  sample 
Interval  of  0.02  second. 

Figure  9.  Same  format  as  Figure  8  for  the  radial  and  vertical  compo¬ 
nents  of  an  explosion  source.  The  bar  Is  1.5E-23  cm  for  a  source  of 
unit  moment. 


C73  SEISMIC  DISCRIMINATION  USING  ONLY  PN 


William  A.  Pepptn 


SUMMARY.  Using  wideband  (0.04  to  10  Hz)  data  recorded  at  the  four  Lawrenre  Liv¬ 
ermore  National  Laboratory  stations,  two  small  Nevada  Test  Site  explosions  (12 
May  1976  1950  GCT;  01  Nov  1978  1525  GCT)  provide  Pn  data  for  this  study.  Using 
vertical  and  radial  components,  the  seismograms  were  Inverted  for  the  seismic 
moment  tensor.  For  each  event,  the  moment  tensor  elements  have  significant  di¬ 
agonal  components,  giving  correct  Identity  of  them  as  explosions.  Moreover,  a 
’master-event"  procedure  Indicates  thct  If  one  event  Is  presumed  to  be  an  explo¬ 
sion,  then  the  other  event  Is  similarly  Identified  with  high  confidence. 

Introduction 

In  view  of  recent  movement  toward  a  test-ban  treaty  Involving  seismic  moni¬ 
toring  at  regional  distances,  It  Is  timely  to  evaluate  new  methods  for  seismic 
discrimination  between  explosions  and  earthquakes  using  regional  seismograms. 
The  recently-developed  seismogram-inversion  method  of  Stump  &  Johnson  (1977)  Is 
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APPENDIX.  Subroutine  to  compute  source  and  receiver  factors,  together 
with  spatial  derivatives  of  the  Green's  functions. 

Variable  Names: 


i 


A 

B 

E 

CL 

ETAA 

ETAB 

ETA-S 

FOR(J) 

GIJK 

LAYER 

OMEGA 

RALEY 


rotated  horizontal  slowness 
.rotated  horizontal  slowness 


"a" 


e,  where 
combined 
vertical 
vertical 
vertical 


b  =  1e  below  (22) 
horizontal  slowness 
slowness  n0  in 
slowness  ng  in 
slowness  at  source. 


in  (17) 
in  (17) 


uwne*: 

n  (6) 
n  (6) 


“i"  in  a  of  (5) 


is  o  or  B 


source  factors  for  point  force  in  direction  Xj 
Green's  function  derivatives  with  respect  to 
layer  number  containing  the  source 
n  defined  below  (7) 

Rayleigh's  function  defined  in  (15) 


X'.:  (see  Johnson,  1974) 


_  _  ^ICiyil  9  lUIIUblUII  UCI  IMCU  Ml 

RECEIV(J)  receiver  factors  for  radial,  transverse,  and  vertical  components 
at  the  free  surface 


i 


i 


/ 


This  subroutine  returns  the  14  non-zero  Green's  function  derivatives 
g..  .  which  are  combined  via  another  subprogram  into  the  canonical  sources. 

Studying  this  subroutine,  the  generalization  to  other  problems  (e.g., 
higher-order  sources;  solutions  within  the  medium)  is  straightforward. 


c 

c 

c 

c 

c 

c 

c 


SUB R 0 IJ  r  I N E  e 3  ( I .  S 0  v  C A  *  E  r  R F L AG  ,  SFLAG  ) 

COMPLEX  LSQfCA 

COMPUTE  GIJfK  AND  RETURN  I  HE  .1.4  NON-ZERO  QUANTITIES 
THROUGH  COMMON  BLOCK  /GREEN/,  THIS  COMMON  BLOCK  IS 
SHARED  BY  FUNCTION  *03*  AND  SUBROUTINE  *R0M84*. 
COMMON  BLOCK  /SUBS/  IS  SHARED  WITH  C0NT02 »  DAPr?,F2, 
ROMCLC r  RTC0E2 r  XFER . 


COMMON  /GREEN/  ETAA  ETA  AS  ,  ETAB  ,  ETABS » GO  ALL , G1 1 1 ,  G13.1. , 

+  G22 1.  ,  G3 1. 1 , G33 1 , G .1. 22 ,  G21 2 , 0232  ,  0322 »  G 1 1 3  r  G 1 33 ,  G223 , 031 3  *  G333 
COMPLEX  ETAA  ,  ET  AAS ,  FT  A1  *  F  T  APS  ,  01. 1  1  ,0131,0271  *  031. 1  ,  0331. « 
f  G  172,021 2, G 2 32?  G 3 2 2  »  C 1 1 3  ,  G 1. 3 3  ••  G  2 2 3  ,  G  3 1 3  ?  0 3 3 •  i 
COMMON  /SUBS/  ALPHA <  7 ) r BETA <  7  >  ,  CMTORS ,  COEF FS ( 20 ) , COSTH » FLAG , 

•}■  INTYPE  f  J  TYPE/ 20)  » NR  0  0  T  ,  K  R  T  ?  L.  V  L  A  Y  E  R  ,  L I S  T  F  L  t ! ..  R  0  M  B  , 
f  NFWC  <  30 )  ,  NEWH  ( 30  )  ,  NEWN  t  NEWNl.  r  N  JTYPE ,  NL.ISI  (20)  >PROD,PU<  /)  ? 

+  r  .  RHO  (.7)  t  ROOT  ,  SHEAR  1  ,  S  I  NTH  ,  T » TCI  HO  R  ,  TFFRMT » TOR  SMC  ,  TS  THOR  ,  L , 
f  UROMB ,  USQ  » XLSQ 
REAL  l  * LROMB * NEWC » NEWH 
COM  PI.  EX  COEFES  ,  PROD  ,  ROOT 
INTEGER  DERI VI 14) , FORCE ( 14 ) ,  COMP  1 1.4) 

COMPLEX  B  »  Cl...  t  DER  ( 3 )  *  V  ACTOR  r  EAC2  ,  FOR  (3)  *  OMEGA  »  RALEY  *  RECE 1 0  <  3 )  , 

■f  T  F  M  F*  ( 1 4  ) 

FQUIVAl  ENCE  (TEMPI  1 )  ,01.11  )  »  (TEMP<2)  >0131  )  »  (TEMP(3>  >0221 )  » 
i  (TEMPI  4) > G31 1 >  » (TEMPIb) ,0331 ) , I  TEMPI A) ,0122) , I » LMP I  7) >0212) > 

+  ( TEMPI  8) »G232) »  I  TEMP  1 9 ) ,8322) * (TEMPI  10) *0113) * (TEMPI  11 ) ,0133) > 
+  (TEMPI  12) ,0223) r (TEMPI  13) ,G313> , I  TEMP  1 14 ) , G333 > 

DATA  DERI V , FORCE , COMP/ 1 ,1,1, 1,1, 2  , 2,2, 3, 3 , 3 , 3,3, 1 ,3,2, 1 ,3, 

+•  2,1  ,3, 2, 1,3, 2,  1,3, 1,1, 2, 3, 3, 1,2, 2 ,3, 1,1, 2, 3, 3/ 


3 


0 

S 

s 


o 

-A 

A 

?>  1 

-  £ 

n  ’"i 

w  3 


CO 


9  C 

as 

HI 


n  n  n  non 


CL=CSQRT  <  LSQ ) 

B=<0.»1* )*E 
DO  30  J— 1 r 1 4 
30  TEMP ( J  >  =0 • 

C 

C  COMPUTE  SOURCE  FACTORS 
C 

IF ( SFLAG .NE . 1 ♦ >  GO  TO  1 
C 

C  P  UP  FROM  SOURCE 
C 

FACT0R=1 .  /CSQRT  ( 2  ♦  *RHO  ( LAYER  )  *ETAAS  ) 

F0R(1>=-CA#FACT0R 

FOR  (2)  =-.B#FACTOR  . 

FOR ( 3 ) =-ETAAS#F ACTOR 
GO  TO  5 

1  IFCSFLAG.NE.3. >G0  TO  2 
C 

C  P  DOWN  FROM  SOURCE 
C 

FACTORS . /CSQRT ( 2 . *RHO ( LAYER ) *ETAAS ) 

FOR ( 1 ) =-CA#FACTOR 

FOR ( 2 ) =-B*F ACTOR 

FOR  <  3 ) =ETAAS#FACTOR 

GO  TO  5 

2  IF<SFLAG.NE.2. )G0  TO  3 
C 

C  SV  UP  FROM  SOURCE 
C 

FACTORS . /CSQRT ( 2 . *RHO ( LAYER  >  *ETABS ) 
FAC2=FACT0R*ETABS/CL 
FOR ( 1 ) =- ( 0 . » 1 . ) *CA*FAC2 
FOR  <  2  >  =- <  0 .  r 1 . >*B*FAC2 
FOR  <  3 )  =  <  0 . r 1 . ) *CL*FAC  TOR 
GO  TO  5 

3  IF<SFLAG.NE.4. )G0  TO  4 

SV  DOWN  FROM  SOURCE 

FACTORS  .  /CSQRT  ( 2 .  *RHO (LAYER  ) #ETABS ) 
FAC2-FAC TGR*E TABS/CL 

FOR ( 1 )  =  <  0 . f 1 . ) *CA*FAC2 
FOR ( 2 )  =  ( 0 ♦  *1.  >*B*FAC2 
FOR <  3 ) -<  0 . » 1 ♦ ) #CL#F ACTOR 
GO  TO  5 

4  IF ( SFLAG. LT .4*5)00  TO  5 

SH  AWAY  FROM  SOURCE 

FACT  0  R - C  UK  C S Q R  T ( 2 . *PU ( LAYER ) *ETABS ) 
FOR ( 1 ) =-B/FACTOR 
FOR ( 2 ) ~C A/ FACTOR 

5  CONTINUE 


-  -  s  :  hj  .  i  iv:  nm 

1  /Ui  •  IS.  -J  TO 
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C 

,  C  RECEIVER  FACTORS*  ONE  FOR  EACH  VALUE  OF  *1*. 

C 

IFIRFLAG.NE.l. >GO  TO  6 

C 

C  P  UP  AT  RECEIVER 
C 

OMEGA=BETA ( X > -2*LSG 

RALEY=0MEGA*0MEGA+4 . *LSG*ETAA*ETAB 

FACT0R=1 . / ( SHEAR1 *RALE Y#CSGRI ( 2 . *RHO ( 1 ) #ETAA ) ) 

RECEIV! 1 )=-FACT0R*4.*ETAA*ETAB#CA 
RECEIV<2)=RECEIV(1)*B/CA 
RECEIV  ( 3 )  =-FACT0R*2  *  *OMEGA*ETAA 
GO  TO  8 

6  IF (RFLAG.NE.2. >G0  TO  7 
C 

C  SV  UP  AT  RECEIVER 
C 

OMEGA=BETA ( 1 ) -2 . *LSQ 

RALEY=0MEGA*0MEGA+4 . *LSG*ETAA*ETAB 

FACT0R=1 . / ( SHEAR 1#R ALE Y#CSGRT < 2 . *RHO < 1 ) *ETAB ) ) 

RECEIV! 1)=FACT0R*<0. *2. > *OMEGA*ETAB*CA/CL 
RECEIV<2)=RECEIV(1 )*B/CA 
RECEIV!3)=-FACT0R*<0. * 4 . ) *CL*ETAA*ETAB 
GO  TO  8 

7  IF(RFLAG.NE.5. )G0  TO  8 
C 

C  SH  UP  AT  RECEIVER 
C 

FACT0R=2 .  /  ( CL.*CSQRT  ( 2 .  *PU  ( 1 )  *ETAB )  ) 

RECE I V ( 1 ) =-B*FACTOR 
RECEIV ( 2 >=CA*FACTOR 

8  CONTINUE 
C 

C  SET  UP  FOR  CORRECT  SPATIAL  DERIVATIVES 
C  "DER*  CONTAINS  FACTORS  FOR  DERIVATIVES  WRT  SOURCE 
C  COORDS  X.YrAND  H.  "SIGNUM"  SET  NEGATIVE  FOR  UPCOMING 
C  WAVES  AND  THE  "H'  DERIVATIVE. 

C 

DERf 1)=CA 
DERI  2) -B 

DER(3)=ETAAS 

IF <  (SFLAG.NE.  1  .  )  ♦  AND .  <  SFLAG . NE  .  3 .  >  ) DL.R <  3 >  =ETABS 

DO  20  J=1 *  14 

SIGNUM=1. 

IF  <  <  < SFLAG , EG . 3 . ) . OR . ( SFLAG . EG . 4 . ) . OR . ( SFLAG . EG . 6 . ) ) . AND . 

+  (DERIV(J) .EG. 3) )SIGNUM-~1. 

IF  <  SFLAG ♦ LT .4.5) GO  TO  21 

IF( (J.NE.l) .AND. (J.NE.3) .AND. < J.NE.&) . AND. ( J.NE.7) .AND. 

+  < J.NE.10) .AND. ( J.NE. 12) )G0  TO  20 
21  TEMP ( J ) =RECEI V ( COMP ( J ) > *FOR < FORCE <  .1 ) >*DER<DERIV< J) )*SIGNUM 
20  CONTINUE 

C  FACTORS  GIJK  RETURNED  THRU  COMMON  ( EQUIVALENCE  TO  TEMP) 

RETURN 

END 


.  .v)*  oa  •  T  to  urs  — ^ 
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ideal  in  two  ways  for  this  problem:  first,  the  method  produces  a  complete  set 

of  source  time  functions,  one  for  each  component  of  the  seismic  moment  tensor 

Mij(t),  and  second,  confidence  analysis  is  Introduced  naturally  in  the  calcula- 
< — > 

tion.  .ne  decision  as  to  the  identity  of  an  event  involves  examination  of  the 
trace  Mil  +  M22  +  M33.  If  this  quantity  is  large  relative  to  the  off-diagonal 

w — *  - - 1  - - ' 

components,  the  event  is  identified  as  an  explosion.  "Large"  can  be  defined  to 
suit  the  criteria  of  the  analysis. 


I  describe  an  effort  to  Invert  the  seismic  phase  Pn  for  the  moment  tensor. 
Two  explosions  are  considered:  12  May  i976,  1950  GCT,  ML  =  4.6,  and  01  Nov 
1978,  1525  GCT,  ML  =  4.3.  The  events  are  selected  because  they  are  fairly  small 
and  because  appropriate  seismic  data  is  available.  Seismic  data  analyzed  con¬ 
sists  of  the  radial,  transverse,  and  vertical  components  of  ground  motion  re¬ 
corded  by  the  wideband  Lawrence  Livermore  National  Laboratory  (LLNL)  stations 
MNV,  200  km  NW,  ELK,,  400  km  NNE ,  KAN,  350  km  ESE ,  and  LAN,  350  km  S  of  the  test 
site.  Several  cons  1  derat i ons  motivate  the  use  of  Pn  alone  in  the  analysis:  (l) 
as  it  is  a  first  arrival,  it  should  be  a  relatively  simple  phase  to  synthesize 
at  regional  distances  as  compared  with,  say  the  complex  Pg  phase;  (2)  it  sam¬ 
ples  a  small  cone  of  Incidence  angles  (l.e.,  near  a  phase  velocity  of  8  kn/s), 
so  provides  a  somewhat  standardized  sampling  of  the  source;  (3)  It  Is  seen  at 
great  distances  in  Asia  where  discrimination  is  of  interest;  (4)  at  low-noise 
sites,  magnitude  4  events  can  probably  be  Investigated  out  to  1,000  km  using 
this  procedure;  and  (5)  as  it  Is  mainly  a  deep-crustal  phase,  we  can  suppose  It 
escapes  some  distortions  caused  by  near-surface  lateral  variations.  The  disad¬ 
vantages  are  that  Pn  is  a  small-amplitude  phase  and  that  it  is  not  a  direct 
wave.  As  a  refracted  wave  it  is  sensitive  not  only  to  first-order  velocity  var¬ 
iations,  but  to  second-order  ones.  In  spite  of  these  difficulties,  the  results 
obtained  appear  quite  promising. 


Method  of  Analysis 


Data  for  the  two  explosions,  recorded  on  analog  magnetic  tape  (FM  mode  de- 
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vlated  40%)  was  digitized  at  60  samplas/s  after  application  of  an  anti-aliasing 
filter  at  10  Hz.  Windows  2.64  seconds  long.  Including  only  the  Pn  and  Its  re¬ 
lated  phases,  but  excluding  direct  arrivals  such  as  PMP,  were  selected  at  the 
four  LLNL  stations.  For  12  May  the  data  traces  used  were*  radial,  transverse 
and  vertical  (R,  T,  Z)  at  MNV;  Z  at  KAN;  R,  T,  and  Z  at  LAN;  Z  at  ELK.  For 
01  Nov  the  same  traces  excepting  the  ELK  vertical  were  used.  The  data  was  con¬ 
ditioned  using  a  10%  cosine  taper.  In  order  to  apply  Stunp  &  Johnson's  (1977) 
procedure,  we  must:  Cl)  select  a  suitable  earth  model,  (2)  compute  the  necessa¬ 
ry  synthetic  seismograms,  and  (3)  Invert  the  seismograms  for  the  moment  tensor 
time  traces  MlJ(t)  using  the  "frequency-domain*  approach.  The  model  selected 
was  a  three-layer  crust  over  a  halfspace  of  7.8  km/s  (Priestley  &  Brune,  1978). 
The  synthetic  seismograms  computed  are  for  a  shallow  (1/2  km)  source.  Ten  com¬ 
ponents  of  motion  must  be  computed,  namely  the  ten  "canonical  source  components" 
(strlkesllp,  dlpsllp,  compensated  linear  vector  dipole,  and  explosion).  Any 
point  source  representable  as  a  linear  combination  of  Green's  function  deriva¬ 
tives  can  be  represented  as  a  linear  combination  of  the  ten  canonical  source 
components  through  the  seismic  moment  tensor  of  order  2  (Stump,  1979).  These 
source  components  were  computed  by  exact  Cagnlard-deHoop  theory  (Peppln,  1981). 
Finally,  the  seismograms  and  source  components  are  Inverted  for  six  time  traces 
MlJ(t),  each  2.54  seconds  long,  using  a  generalized  Inverse  technique. 


Results 


In  Figure  1  we  compare  observed  vertical-component  Pn  waveforms  with  the 

synthetic  seismogram  for  an  explosion  source  with  a  stepup  time  of  0.3  second. 

As  can  be  seen,  the  data  are  distressingly  unlike  the  synthetics.  Thus,  the 
outlook  for  success  appears  quite  bleak.  We  fit  only  the  vertical  and  radial 
components,  downweighting  the  transverse  ones,  because  the  transverse  components 
are  hopelessly  contaminated  by  multlpathlng  (synthetic  transverse  two  orders  of 
magnitude  smaller  than  verticals,  but  observed  transverse  of  comparable  ampli¬ 
tude).  Examples  of  the  fits  between  observed  and  synthetic  seismograms  are 
given  In  Figure  2.  It  Is  not  notorious  that  the  fits  are  very  close,  because 
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the  nuiiber  of  components  being  fit  Is  small  (seven  for  01  Nov  and  eight  for  12 
May)»  given  as  few  as  six  components,  an  exact  fit  would  be  obtained  for  almost 
any  Green's  function  selected.  Note  In  the  figure  the  tremendous  difference  In 
signal  character  for  the  two  explosions;  the  Inversion  simply  places  this  char¬ 
acter  In  the  time  traces  of  MlJ(t),  giving  an  unacceptably  long  and  complex 
source  time  history  for  these  events. 


Consider  now  the  maximum  (unsigned)  trace  amplitudes  of  MlJ(t)  for  the  two 
events.  They  are  (normalized  to  unit  determinant) 
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.  will 

for  the  12  May  event  and 
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for  the  01  Nov  event.  We  see  that  the  diagonal  elements  dominate,  which  sug¬ 
gests  that  each  event  should  be  Identified  as  an  explosion.  Moreover,  suppose 

we  declare  that  the  12  May  event  is  Indeed  an  explosion.  Let  factors  P1J  be  de¬ 
termined  such  that 

S 

M1J  (12  May)  ♦  Pi J  -  de^lj. 


X  ui 


and  let  P1J  be  added  to  the  MiJ  array  for  the  01  Nov  event. 

j  —j 

01  Nov  array  becomes  almost  pure  explosion  In  character: 

— 
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Then  the  adjusted 
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Thus,  a  ’master-event *  method  Is  able  to  establish  quite  clearly  the  Identity  of 
an  unknown  If  a  known  Is  established.  This  Is  a  significant  result  because: 

<1)  the  signal  character  of  the  two  events  differs  so  greatly  (Figure  2>;  (2) 

12  May  Included  data  at  Elko  while  01  Nov  did  not;  (3)  12  May  was  In  Yucca  Val¬ 
ley  and  01  Nov  was  on  Pahute  Mesa,  and  (4)  01  Nov  was  overburled. 

It  Is  of  some  Interest  to  Investigate  the  devlatoric  part  of  the  Mlj(t).  L_i 
At  any  speclfiec  time  this  array  (with  signs)  has  a  set  of  eigenvectors  which 
point  In  the  principal  directions  of  stress  (Gilbert,  1973;  Stump,  1979). 

Consider  a  streographlc  plot  of  the  eigenvectors  for  the  12  May  event  associated 
with  the  minimum  and  maximum  eigenvalues  (axes  of  compression  and  extension), 
shown  In  Figure  3.  These  axes  cluster  mainly  about  two  areas  on  the  stereonet. 
Strikesllp  motion  Is  Indicated  In  a  sense  similar  to  some  of  the  mechanisms  ob¬ 
served  by  Hamilton  &  Healy  (1969)  for  BENHAM  aftershocks  on  Pahute  Mesa. 

Di scussi on 

The  results  agree  In  one  Important  aspect  with  other  efforts  at  Inversion 
(Stump,  1979;  Stump  &  Johnson,  1981).  It  appears  that  the  amplitudes  of  the 
Mlj  will  be  Insensitive  to  errors  in  the  Green's  functions  selected,  while  the  l— J 

phase  information  is  Indeed  quite  sensitive.  Thus,  we  can  place  little  stock  In 
the  extended  duration  of  the  source  (no  doubt  an  artifact  of  Incomplete  knowl¬ 
edge  of  the  Pn  phase);  nor  can  we  place  much  confidence  In  the  Inferred  direc¬ 
tions  of  stress  for  the  devlatoric  part  of  the  solution.  However,  It  Is  clear 
that  with  better  knowledge  of  the  structure  these  results  could  be  Improved. 
Furthermore,  the  results  presented  here  based  only  on  amplitude  of  the  Mlj  ap-  l— J 

pear  to  be  of  considerable  interest. 

Note  that  we  have  attempted  no  Inversion  of  earthquake  data  using  the  Pn 
phase  (lack  of  any  known  suitable  event  near  the  test  site  recorded  on  the 
three-component  LLNL  stations).  We  have  good  reason  to  suppose  that  any  Pn  data 
will  Invert  to  a  significant  Isotropic  component  In  Mlj,  since,  for  example,  Pn  u 
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contains  no  significant  SH  component  (apart  from  multipathlng)  which  would  tend 
to  drive  down  the  Isotropic  component.  However,  It  Is  difficult  to  believe  that 
our  master-event  technique,  applied  to  an  earthquake,  would  give  so  nearly  an 
Isotropic  result.  Moreover,  In  Figure  3  we  see  the  potential  for  analysis  of 
the  devlatorlc  part  of  the  source.  For  a  pure  explosion,  we  should  have  expect¬ 
ed  a  patternless  distribution  of  eigenvectors  (this  Is  the  degenerate  eigenvalue 
problem).  Figure  3  provides  the  opportunity  to  estimate  the  focal  mechanism 
from  a  considerable  number  (up  to  one  estimate  for  each  seismogram  time  sample, 
here  128)  of  independent  estimates.  In  sum,  this  line  of  analysis  appears  de¬ 
finitely  worth  pursuing  in  connection  with  seismic  discrimination,  with  the  pre¬ 
sent  results  quite  encouraging. 
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Figure  Captions 

Figure  1.  Comparison  of  synthetic  selsmograsm,  vertical  component,  with 
observed  vertical  components  at  Landers  (left)  and  Elko  (right).  Top 
traces?  computed  ground  displacement,  step-function  source;  second 
traces!  after  applying  Instrument  correction,  essentially  to  ground  velo¬ 
city;  third  traces!  convolution  of  second  trace  with  a  source  time  func¬ 
tion  of  0.3  second  stepup  time;  bottom  traces!  data.  The  synthetics  con¬ 
sist  almost  entirely  of  the  downgoing  Pn  and  its  reflection  at  the  free 
surface;  other  phases  are  multiples  1n>the  top  (3  km  thick)  layer. 


Figure  2.  Comparison  of  observed  seismograms  (top  of  each  pair)  with  fit¬ 
ted  seismograms  (bottom  of  each  pair).  Left!  12  May  1976;  right!  01  Nov 
1978.  On  the  radial  traces  away  from  the  source  is  up,  and  on  the  vertical 
traces  compression  Is  downward.  The  relative  scales  of  the  pairs  are  cor¬ 
rect.  The  fits  are  as  good  on  vertical  and  radial  traces  not  shown,  but 
not  on  the  transverse,  which  were  downweighted. 


Figure  3.  Stereo  projection  of  axes  of  compression  (left)  and  extension 
(right).  Each  point  gives  the  eigenvectors  of  Mlj(t>  at  a  particular  time 
step.  The  inferred  mechanism  is  NW  strikeslip,  either  right-lateral  or 
lef t-lateral ,  depending  on  how  the  clusters  on  the  stereonet  are  selected. 
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